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1 . 0 INTRODUCTION 


1 . 1 Background 

Sensitivity  analysis  is  receiving  a widening  interest  because  of  its  role 
in  model  validation  and  the  need  to  evaluate  model  results  in  terms  of 
variations  of  input. 

The  complexity  and  sophistication  of  contemporary  models  and  the 
attendant  justification  of  model  results  in  matters  of  policy  and  planning 
have  created  a demand  for  a wide  range  of  sensitivity  methodology. 

Classical  sensitivity  analysis,  which  was  formalized  in  the  last  century, 
is  no  longer  adequate  to  meet  the  requirements  of  most  modelling  tasks.  It 
was  based  either  on  an  application  of  statistical  methodology  to  observed  data 
or  on  the  use  of  differentials  of  variables  of  interest.  These  approaches  are 
still  used,  particularly  in  experiments  of  an  engineering  nature  or  in  the 
context  of  a controlled  laboratory  protocol,  but  they  do  not  meet  the  needs, 
for  various  reasons  which  will  be  discussed,  of  contemporary  model  designs 
intended  to  be  run  on  high-speed  computers.  The  technical  literature  is 
reporting  new  analytical  and  statistical  methods  for  validation  and 
sensitivity  analysis.  These  methods  employ  techniques  which  were  developed 
specifically  for  the  complex  models  of  policy  and  allocation  problems. 

Sensitivity  analysis  is  an  important  subject  in  model  evaluation  because 
it  (1)  aids  in  validating  a model  by  contributing  to  an  understanding  of  the 
correspondence  between  a model  and  its  subject,  (2)  assists  in  the 
identification  of  variables  which  contribute  only  marginally  to  functional 
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relationships  among  variables  or  elements  in  the  model,  (3)  measures  the 
quantitative  impact  of  variations  of  input  on  model  results,  (4)  tests  the 
sensitivity  of  the  level  of  representation  of  activities  in  the  model  to 
variations  which  effect  actual  perturbations  or  changes,  and  (5)  contributes 
to  measuring  the  degree  to  which  input  data  may  be  altered  without 
significantly  exceeding  a specified  range-of-interest  in  the  output. 

A sensitivity  analysis  of  a computational  process  or  model  may  have  to 
deal  with:  (1)  errors  which  are  inherent;  that  is,  they  are  assumed  to  exist 
because  of  sampling  or  collection  protocols  (2)  errors  which  vary 
stochastically;  and  (3)  errors  which  are  imposed  upon  the  process  by  external 
circumstances  or  controls.  The  first  kind  of  error,  those  which  are  embedded 
in  the  data,  contribute  to  errors  in  the  output  primarily  as  displacements  or 
scaling  factors.  The  second  kind  of  error,  those  which  are  random,  project  or 
convolute  their  randomness  on  the  output,  thereby  giving  the  output  the 
character  of  a statistical  distribution.  The  third  kind  of  error,  that  which 
is  imposed,  creates  a displacement  error  in  the  output.  The  imposed  error  is 
evaluated  by  methods  used  to  ascertain  system  biases.  The  advanced 
methodology  of  sensitivity  analysis  for  evaluating  these  error  types  are  drawn 
from  specialized  applications  of  mathematics,  physics,  and  statistical 
sampling. 

1.2  Motivation  and  Objectives 

Sensitivity,  in  the  sense  of  responses  to  variations,  is  of  two  broad 
kinds:  qualitative  and  quantitative. 

Qualitative  sensitivity  is  used  in  this  paper  to  denote  the  capability  of 
a model  to  respond  to  dynamic  changes  in  the  subject  being  modeled.  A closely 
related  definition  is  the  level  of  representation  of  the  subject  by  its  model. 


2 


The  sensitivity  referred  to  here  is  the  difference  between  the  representation 
and  the  subject,  and  the  significance  of  that  difference  in  terms  of  second 
and  third  order  effects. 

An  example  of  dynamic  qualitative  sensitivity  would  be  a model  of  a 
queueing  problem  in  which  the  logic  of  the  model  permits  a dynamic  reponse  to 
queue  sizes  (demand)  by  opening  or  closing  channels  (service)  according  to  the 
circumstances.  Another  example,  though  not  commonly  used  in  typical  model 
design,  is  a heuristic  model  which  alters  its  process  based  on  a "learning" 
ability  or  the  use  of  artificial  intelligence  to  select  or  alter  its  logic 
according  to  accumulated  information.  The  use  of  "dynamic"  should  not  be 
confused  with  an  association  of  time  as  a model  variable,  although  it  is  not 
excluded  from  the  discussion  above  - dynamic,  as  used  here,  refers  to  the 
ability  of  the  model  to  automatically  modify  its  logic  or  computational 
process  according  to  critical  perturbations  in  specified  state  variables  which 
may  or  may  not  be  functions  of  time.  Most  models  deal  with  a specified 
scenario  with  almost  no  capability  to  simulate  potential  system  dynamics.  It 
is,  therefore,  difficult  to  validate  these  models  or  to  understand  thoroughly 
their  range  of  application. 

Quantitative  sensitivity  is  defined  as  the  numerical  measure  of  changes 
in  output  to  variations  of  input.  This  is  the  traditional  view  of  sensitivity 
analysis . 

The  motivation  for  preparing  a survey  of  sensitivity  methodology 
addresses  both  the  need  to  recognize  and  analyze  qualitative  and  quantitative 
sensitivity  problems.  The  qualitative  aspect  addresses  model  design;  the 
quantitative  aspect  addresses  responses  to  parametric  variation. 
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The  motivation  and  objectives  of  this  survey  are  interlinked  in  the  need 
and  desire  to  provide  sensitivity  methodology  and  reference  in  a single 
source;  to  advertise  to  the  NBS  staff,  and  others,  the  importance  of  the 
subject;  to  improve  the  understanding  of  sensitivity  analysis  in  the  design  of 
models;  and  to  enhance  the  ability  to  perform  a more  efficient  model 
validation  through  the  use  of  qualitative  and  quantitative  sensitivity 
methodology. 

This  survey  was  written  to  provide  a general  overview  of  sensitivity 
analysis  to  the  model  developer  arid  is  not  intended  to  provide  a precise 
theoretical  and  comprehensive  treatment  of  the  methods  included  in  the  survey. 
1.3  Organization 

The  paper  is  organized  into  six  sections  which  develop  and  evaluate  the 
classical  and  modern  methods  of  sensitivity  analysis.  Section  2.0,  Error 
Types  and  Model  Structures,  defines  the  error  types  incurred  in  modeling,  and 
in  statistical  and  physical  experiments;  Section  3.0,  Classical  Methods  of 
Sensitivity  Analysis,  presents  the  familiar  theory  of  errors  and  numerical 
significance;  Section  4.0,  Statistical  Methods  of  Sensitivity  Analysis, 
describes  the  basic  statistical  approaches  to  measuring  model  performance  and 
other  aspects  of  statistical  interest;  Section  5.0,  Analytical  Methods  of 
Sensitivity  Analysis,  presents  an  assortment  of  methods  which  constitute  the 
major  advances  in  different  lines  of  sensitivity  methodology;  Section  6.0, 
Special  Analytical  Methods  of  Sensitivity  Analysis  presents  some  very 
interesting  techniques  and  concepts  for  certain  kinds  of  sensitivity  problems 
with  broad  application;  and  Section  7.0,  Modeling  and  Sensitivity  Analysis,  in 
which  methods  of  model  assessment  and  design  concepts  are  presented  and 
evaluated . 


4 


Although  it  may  not  be  possible  to  discuss  every  method  presented  with 
equal  knowledge  or  to  know  the  consequences  of  its  application,  an  attempt  to 
evaluate  some  of  the  methods  has  been  made.  The  evaluation,  when  given,  is 
based  on  information  from  cited  references  which  includes  comments  on 
efficiency,  ease  of  use,  cost,  preparation,  sample  size,  and  output 
description. 

Equations  cited  within  quoted  material  retain  their  source  designations; 
equations  in  the  text  are  numbered  according  to  the  major  section  in  which 
they  appear. 
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2.0  ERROR  TYPES  AND  MODEL  STRUCTURES 


2.1  Introduction 

This  section  introduces  the  basic  types  of  errors  incurred  in  data  and 
model  analysis.  It  presents  the  context  of  these  error  types  by  discussing 
generic  functional  forms  on  which  sensitivity  analysis  is  carried  out,  and 
rounds  out  the  discussion  of  these  subjects  with  a description  of  model  types 
and  their  relation  to  sensitivity  analysis.  Taken  together  these  three 
subjects  define  fundamental  relationships  among  the  model,  the  representation 
of  the  subject,  and  the  data.  It  is  from  these  relationships  that  data  and 
model  uncertainties  have  their  origins. 

2.2  Error  Types 

There  are  five  major  types  of  errors  ranging  from  the  accuracy  of  input 
data  to  the  implicit  errors  introduced  by  the  selection  of  computational 
methods,  algorithms,  and  the  formulas  themselves  which  are  used  to  evaluate 
model  sensitivity. 

2.2.1  Accuracy  of  Numerical  Values 

Although  there  may  be  many  ways  in  which  uncertainties  can  be 
incorporated  into  input  data,  the  error  which  is  emphasized  here  is  that  of 
the  number  of  significant  digits  in  data  of  a cohort  or  homogeneous  file. 
Significance  of  numbers  is  a field  of  numerical  analysis  which  has  been 
studied  for  a long  time  and  most  of  the  rules  addressing  numerical  precision 
are  well  known.  Several  of  these  rules  will  be  summarized  here  to  refresh  our 
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recognition  of  the  subject  and  to  refer  to  these  rules  when  our  examination 
reaches  the  study  of  the  propagation  of  errors  in  a computatonal  process. 
Following  Reference  [1],  Kelly  gives  three  general  rules  dealing  with 
numerical  precision: 

(1)  If  the  first  significant  figure  of  a number  is  k,  and  the  number  is 
correct  to  n significant  figures,  then  the  relative  error  in  this  number  is 
less  than 

1 (2.1) 

klOn_1 

where  the  relative  error  is  defined  as  the  error  divided  by  the  number. 

(2)  If  the  relative  error  in  a number  is  less  than 

1 

(k+1 ) 10n“l  (2.2) 

the  number  is  correct  to  n significant  figures  or  is  in  error  by  less  than  a 
unit  in  the  nth  significant  figure. 

(3)  A number  is  correct  to  n significant  figures  if  the  relative  error  of 
the  number  is  not  greater  than 

1 1 (2.3) 

2 • 10n 

The  propagation  of  errors  in  significant  digits  depends  entirely  on  the 
distribution  of  the  errors  throughout  the  data  set  and  the  computational 
process  to  which  the  data  are  subjected  in  the  model. 

2.2.2  Process  Errors 

Process  errors  are  those  errors  introduced  by  rounding  and  by  truncation. 
These  errors  arise  primarily  from  register  lengths  in  computer  hardware,  and 
are  entirely  dependent  on  the  mechanical  features  of  the  equipment  used  for 
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the  model  computation.  They  may  or  may  not  be  important  to  sensitivity 
analysis  since  most  computers  provide  floating  point  calculations.  It  is  not 
unreasonable  to  make  the  point  that,  for  most  simulations,  the  input  data  are 
far  less  accurate  than  the  capacity  for  accuracy  provided  by  most  modern 
computers.  However,  round-off  error  is  an  extremely  important  subject  and  has 
been  treated  widely  in  the  literature  as  it  relates  to  and  impacts  on  the 
accuracy  of  repeated  calculations,  convergence  criteria,  and  decisions  based 
on  numerical  differences. 

2.2.3  System  Errors 

There  are  several  categories  of  errors  which  can  be  thought  of  as 
associated  with  system  processes:  those  that  come  from  system  bias,  those 

that  arise  from  system  noise,  those  that  originate  from  measurement,  those 
that  occur  in  data  generation  protocols,  and  those  that  arise  from  data 
sampling.  These  may  be  regrouped  under  three  more  convenient  categories: 
systematic  errors;  to  cover  system  bias  and  noise;  measurement  errors , to 
cover  the  errors  of  measurement,  sampling,  and  generation;  and  structural 
errors , to  describe  errors  imposed  by  the  computational  process  of  the  model. 
It  is  necessary  to  examine  the  possibility  that  the  structure  of  the  model  or 
experimental  design  may  directly  affect  the  results  of  an  experiment.  This 
point  involves  both  qualitative  sensitivity  and  quantitative  sensitivity;  it 
is  an  important  aspect  of  sensitivity  analysis  because  model  design  may  impose 
an  error  of  perception  or  bias  on  the  structure  of  the  model.  The  general 
theory  of  system  errors,  as  they  occur  in  the  context  of  physical  experiments, 
is  nicely  covered  in  References  [3],  [4],  and  [5]. 
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2.2.4  Computational  Errors 

In  this  class  of  errors  there  are  three  levels  at  which  qualitative  or 
quantitative  errors  are  likely  to  be  incurred.  A qualitative  error  is  an 
error  imposed  by  selecting  a particular  method  of  computation  or 
representation;  and  a quantitative  error  is  an  error  associated  with  numerical 
calculations . 

The  first  level  the  source  of  error  is  the  method  of  computation, 
representation,  or  degree  of  aggregation  of  the  subject.  By  choosing  a 
particular  process  for  a computation,  approximations  to  the  subject  are 
automatically  imposed  which  invest  the  process  with  attendant  quantitative 
errors.  The  choice  of  algorithm  or  model  logic  is  a qualitative  decision,  but 
when  the  choice  is  made,  quantitative  errors  are  imposed  on  the  results 
obtained  from  the  method  of  computation. 

The  second  level  of  error  is  that  which  originates  from  the  equations 
used  by  the  selected  method  of  computation.  This  level  of  error  is 
characterized  by  the  calculation  of  quantitative  errors  in  the  form  of 
relative  functional  errors,  sensitivity  coefficients  (discussed  in  detail  in 
3. 3. 1.4),  and  parametric  variations. 

The  third  level  of  error  is  the  error  of  calculation,  which  results 
from  a sequence  of  arithmetic  operations  upon  numbers  in  the  equations  of  the 
method  of  computation.  This  level  addresses  issues  of  significance  of  the 
input  data,  significance  of  results  of  arithmetic  operations,  and  the  net 
precision  of  the  final  result  of  the  computational  process.  Each  level 
contributes  to  the  propagation  of  errors  in  either  the  quantitative  or 
qualitative  sense,  or  both. 
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2.2.5  Errors  of  Sensitivity  Analysis 


The  methods  of  sensitivity  analysis,  whether  they  be  statistical  or 
variational  formulations,  or  any  of  the  many  forms  described  in  this  report, 
may  themselves  create  errors  because  some  are  approximations  or,  by  the 
argument  of  Section  2.2.4  above,  they  may  induce  errors  or  biases  because  they 
involve  sequences  of  calculations. 

The  methods  of  quantitative  sensitivity  analysis  are  primarily  either 
statistical  in  nature,  use  local  approximations  (Taylor  series),  or  employ  a 
closed  functional  solution  to  obtain  the  desired  sensitivity  coefficients. 
There  also  may  be  assumptions  which  bear  on  numerical  characteristics,  such  as 
distribution  or  error  magnitude,  which  contribute  to  the  total 
quantitative  sensitivity  of  a model's  computational  process. 

2.3  Function  Types 

Error  analysis  is  ultimately  contingent  on  the  nucleus  of  arithmetic 
operations  embedded  in  a method,  a model,  or  an  experiment.  These  operations 
collect  and  accumulate  qualitative  and  quantitative  errors  which  are  the 
pursued  objects  of  sensitivity  analysis.  A model  is  often,  but 
simplistically , described  as  a function,  as  in 

y = f( a ,x)  (2.4) 

where  the  "a"  is  a set  of  input  control  data  and  "x"  is  a set  of  input  state 
data.  The  model  is  represented  by  the  operator  ”f"  which  acts  on  "a"  and  "x” 
to  produce  output  "y” . Equation  (2.4)  is  the  basic  representation  which 
is  used  in  sensitivity  analysis  for  the  development  of  absolute  and  relative 
errors  in  the  variable  y. 

An  important  extension  of  this  concept  is  the  system  which  involves 
embedded  functional  relationships,  as  in 

y = f(  a ,x) 

u = g(b,y,v)  (2.5) 


10 


in  which  "b"  is  a set  of  control  input  data,  y is  the  output  from  the  previous 
equation,  "v"  is  a set  of  control  state  variables  which  may  include  some 
subset  of  "x",  and  "g"  is  an  operator.  This  process  of  chaining  is  common  to 
most  simulations  and  models  of  any  degree  of  complexity. 

The  next  step  up  in  functional  formulation  is  the  system  of  simultaneous 
equations,  given  as 

y = Ax  (2.6) 


for  a linear  system,  where  A is  a square  matrix  and  x,  y are  vectors;  or  as 

y’  = f(x,y) , (2.7) 

the  differential  form,  with  x,y  as  vectors  and  y'  as  the  derivative  of  y with 
respect  to  x.  Equation  (2.7)  could  also  be  expressed  as 

y = f(x,y ,t) , (2.8) 

where  y is  the  derivative  of  y with  respect  to  time,  "t”.  Equation 
(2.8)  is  the  general  form  for  the  sensitivity  function,  discussed  below  in 
Section  5.6.  Equation  (2.6)  appears  in  many  problem  formulations  and  most 
particularly  as  the  basic  form,  with  an  objective  function,  of  a linear 
program.  If  in  equation  (2.6)  y = Xx,  then  the  familiar  eigenvalue 
formulation  emerges. 

2. A Model  Types 

It  is  customary  to  refer  to  model  types  as  descriptive , prescriptive , or 
predictive , and  within  these,  the  models  as  being  static  or  dynamic.  These 
categories  designate  purpose  and  structure,  but  they  are  not  particularly 
useful  distinctions  for  the  purposes  of  sensitivity  analysis. 
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As  this  Section  has  attempted  to  indicate,  there  is  a progression  of 
error  analysis  which  spans  the  numerical  significance  of  data  at  the  lowest 
level  to  the  extensive  variety  and  complexity  of  models  and  simulations  at  the 
highest  level.  At  each  stage  of  this  progression  the  complexity  of  the  error 
analysis  increases  in  degree  over  the  previous  stage,  which  requires,  in  turn, 
an  increase  in  the  sophistication  of  the  sensitivity  methodology  as  the 
progression  moves  from  simple  to  complex  model  processes. 

Many  of  the  methods  of  sensitivity  analysis  have  been  developed  to 
evaluate  sophisticated  models  of  complex  subjects,  and  therefore  are,  in  many 
cases,  designed  for  a particular  treatment  or  mathematical  formulation.  There 
are,  then,  at  least  seven  general  types  of  model  processes  that  can  be 
identified;  and  each  has  its  own  particular  method  of  sensitivity  analysis: 

(1)  independent  algebraic  or  differential  equations,  (2)  chained  dependent 
algebraic  equations,  (3)  system  of  algebraic  equations  (4)  system  of 
differential  equations,  (5)  system  of  stochastic  differential  equations,  (6) 
statistical  distribution  functions  or  variables,  and  (7)  the  stochastic  or 
deterministic  iterated,  closed  program  (algorithm).  There  are  certainly  many 
other  variations  of  these  types  that  could  be  listed,  but  to  dwell  on 
categories  would  obscure  the  principal  point,  namely,  that  sensitivity 
methodology  consists  of  general  techniques  for  analyzing  subject  formulations 
rather  than  analyzing  specific  model  types.  Even  the  model  itself,  because  of 
its  organization  of  components,  degree  of  complexity,  and  level  of  aggregation 
may  be  considered  as  a process  which  is  appropriate  for  sensitivity  analysis 
as  a part  of  the  larger  process  of  validating  subject-model  correspondence. 
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Model  types,  therefore,  which  relate  to  sensitivity  analysis  should  be 
categorized  in  terms  of  the  representation  of  the  subject.  This  concept  is 
essentially  the  reverse  of  the  conventional  approach,  which  is  to  identify  the 
model  by  the  purpose  of  the  model  and  the  use  of  its  output. 
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3.0  CLASSICAL  METHODS  OF  SENSITIVITY  ANALYSIS 


3.1  Introduction 

The  principal  subjects  associated  with  classical  sensitivity  analysis  are 
computational  precision  (primarily  rounding  errors  and  numerical  error 
estimation) , the  general  theory  of  errors  (addressing  variations  of  functions 
and  products  or  quotients  of  random  variables) , and  differential  formulations 
and  approximations  (which  provide  closed,  analytical  expressions  of  error  and 
sensitivity).  These  subjects  are  examined  in  the  following  sections.  The 
material  is  not  intended  to  be  exhaustive  but  rather  in  the  spirit  of  a survey 
of  issues  of  primary  interest  to  model  developers  rather  than  to  numerical 
analysts . 

3.2  Precision  In  Computation 

Computational  precision  is  represented  primarily  by  the  net  numerical 
significance  of  arithmetic  calculations  and  by  the  analysis  of  roundoff  errors 
in  computational  operations. 

3.2.1  Numerical  Significance 

Section  2.2.1,  Accuracy  of  Numerical  Values,  introduced  several  rules 
about  relative  errors  in  numbers  and  measures  of  these  errors.  There  are  two 
ways  of  defining  a significant  figure:  (1)  as  any  of  the  ten  digits,  except 

zeros  which  are  used  only  to  locate  the  decimal  point,  and  (2)  any  digit  which 
signifies  the  amount  of  quantity  in  the  place  in  a number  in  which  it  stands. 
We  now  state  from  Kelly,  Reference  [1],  the  relation  of  significance  to  simple 
arithmetic  operations. 

(1)  If  two  numbers  are  added  or  subtracted  the  maximum  error  is  the  sum 
of  the  maximum  errors  of  the  two  numbers, 

(2)  If  two  numbers  are  multiplied  or  divided,  the  result  has  a maximum 
relative  error  equal  to  the  sum  of  the  relative  errors  of  the  two  numbers. 
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In  addition,  experience  has  shown  that  errors  of  a sum  or  multiplication 
compensate  to  a degree,  especially  in  a long  string  of  calculations.  The  loss 
of  significance  by  subtraction  is  the  principal  source  of  error. 

Additional  rules  advise  on  accuracy  where  a variation  is  introduced: 

(3)  retain  enough  significant  figures  in  a number  to  include  the  place  in 
which  the  least  significant  figure  of  the  variation  occurs;  and 

(4)  in  adding  or  subtracting  several  numbers  deal  with  absolute  error, 
not  the  relative  error, 

(5)  when  several  values  are  to  be  multiplied  or  divided  use  the  largest 
relative  error  in  computing  the  number  of  significant  digits  of  the  result. 
Stoer  & Bulirsch,  Reference  [2] , have  written  a comprehensive  and  modern 
treatment  of  computer  arithmetic,  error  propagation  and  interval  analysis, 
which  complements  the  older  texts  on  numerical  analysis. 

3.2.2  Rounding  Errors 

This  subject  is  usually  included  under  discussions  on  computational 
precision  where  the  roundoff  is  made  to  a number  in  a particular  digit  of 
the  number,  rounded  up  or  down  in  that  digit,  and  the  rest  of  the  number  is 
truncated.  There  is,  however,  a more  sophisticated  analysis  of  the  errors 
associated  with  rounding,  which  is  given  in  Reference  [2],  Chapter  1,  and  is 
there  referred  to  as  statistical  roundoff  estimates.  This  concept  defines 
the  relative  roundoff  error  which  results  from  an  elementary  operation  as  a 
random  variable  on  a given  interval  of  values.  The  roundoff  errors  are 
assumed  to  be  independent  with  a given  distribution  of  values.  Seen  in  this 
way,  the  evaluation  of  roundoff  errors  proceeds  along  classical  statistical 
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determinations  of  the  means  and  variances  of  results  of  computations.  The 
assumption  of  independence  is  critical:  if  the  roundoff  errors  are  not 

independent  the  error  analysis  requires  a more  complicated  formulation. 

3.3  General  Formulations  of  Error  Analysis 

This  section  presents  the  subjects  which  comprise  the  bulk  of  classical 
sensitivity  analysis.  The  familiar  results  of  numerical  analysis,  statistical 
theory  and  error  estimates  are  developed  and  evaluated.  These  methods,  and 
the  statistical  methods  described  in  Section  3.3.2  below,  are  the  most  widely 
used  of  all  the  techniques  available  to  model  assessors.  A summary  of  the 
principal  formulas  of  the  propagation  of  errors  and  classical  sensitivity 
analysis  is  provided  in  Appendix  A. 

3.3.1  Theory  of  Errors 

Within  this  broad  subject  there  are  three  subordinate  subjects  which  are 
discussed  separately  in  order  that  their  contributions  to  the  general  theory 
of  errors  are  clearly  delineated.  These  three  subjects  are  product  and 
quotient  analysis;  errors  of  integration,  interpolation  and  differentiation; 
and  propagation  of  errors. 

3. 3. 1.1  Quotient  and  Product  Analysis 

There  are  two  basic  relationships  in  the  theory  of  errors  which  are 
fundamental  to  sensitivity  analysis.  These  relationships  are  the  statistical 
expression  of  the  means  and  standard  deviations  of  the  product  and  quotient  of 
two  random  variables. 
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Given  that  x and  y are  random  and  normally  distributed  with  means  yx  and 
Uy,  and  standard  deviations  ox  and  Oy,  let 

w = x/y  , (3.1) 

then  the  mean  of  w,  w is: 

w = yx/yy  , yy  * 0.  (3.2) 


The  standard  deviation  of  w was  developed  in  a basic  paper  by  Fieller 
[6],  and  is  as  follows:  If  w = x/y,  and  x and  y are  independent  and  normally 
distributed  with  means  yx  and  Py  and  standard  deviations  ax  and  Oy,  then  w 
will  be  approximately  normally  distributed  with  mean  given  in  (3.2)  above  and 
a standard  deviation  of: 


aw  y 


Hx 

y 


’x 


+ fz-  - 2r 


1/2 


(3.3) 


b x2  W y 2 y y 

where  r is  the  correlation  coefficient  and  Py/cfy  > 5.  See  Appendix  A,  Part 
II,  for  more  details  on  this  theory. 

Craig  [7]  developed  a similar  formula  for  the  product  of  random 
variables 

w=xy.  (3.4) 

If  x and  y are  independent  and  normally  distributed  with  means  yx  and  yyi 
and  standard  deviations  ox  and  Oy,  then  the  mean  of  w is 

w = yxpy  , (3.5) 

and  the  standard  deviation  is: 

aw  = (Mx^Oy^+yy^ax^+ox^oy2)l/2  (3.6) 


In  the  case  that  x and  y are  not  independent,  but  normally  distributed, 
°w  is  (from  [7],  page  9): 


® w [ P x^^y^i~y  2ryxPyOxay-l-0x2ay2  ( 1+r2 ) ] 


1/2 


(3.7) 
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where  r is  the  correlation  coefficient.  Reference  [8]  gives  a procedure  for 
obtaining  either  one-sided  or  two-sided  tolerance^  limits  of  the  product  and 
quotient  of  independent  normal  variates.  It  is  assumed  in  this  procedure  that 
the  distributions  have  positive  means  and  small  coefficients  of  variation 
(ox/yx).  Huntington  [9]  derives  frequency  distributions  for  the  product  and 
quotient  of  two  independent  random  variables,  and  Curtiss  [10]  provides  a 
general  and  thorough  development  of  the  frequency  and  distribution  functions 
of  the  quotient  of  two  random  variables.  (Curtiss  concludes  his  paper  with  a 
brief  treatment  of  the  distribution  of  the  product  of  two  random  variables.) 
His  theorems  are  generalized  and  expressed  as  integral  functions  depending  on 
the  positive  and  negative  domains  of  the  product. 

3. 3. 1.2  Integration,  Interpolation,  Differentiation  Errors 

These  mathematical  methods  are  well  known  and  thoroughly  treated  in  texts 
on  numerical  analysis.  The  advent  of  high-speed,  sophisticated  computers  has 
made  it  possible  to  build  models  of  extreme  complexity.  This  capacity  has 
required  a new  development  of  numerical  procedures,  and  as  a consequence  the 
classical  numerical  techniques  have  been  revitalized  or  replaced  by  advanced, 
innovative  methods. 

No  attempt  will  be  made  here  to  restate  the  classical  methodology  nor 
will  any  attempt  be  made  to  reproduce  any  of  the  many  advanced  numerical 
techniques  except  to  recognize  their  existence  and  their  immense  importance  to 
many  diverse  areas  of  application.  Software  packages  abound  which  contain 

^Not  to  be  confused  with  confidence  limits.  Tolerance  limits  are  defined  as 
the  limits  which  include  at  least  "p"  percent  of  the  observations  at  some 
prescribed  level  of  confidence. 
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specific  and  generic  numerical  tools  for  virtually  every  computational  need, 
and  the  underlying  theory  of  these  developments  have  given  numerical  analysis 
an  enhanced  and  extensive  revitalization. 

3. 3. 1.3  Propagation  of  Errors 

The  theory  of  propagation  of  errors  should  be  more  properly  called  the 
theory  of  small  errors.  It  was  developed  about  a century  ago  and  most  of  the 
early  texts  which  appeared  around  the  turn  of  the  century  dealt  with  analysis 
of  observations  generated  from  physical  experiments  and  calibration  of 
instruments.  More  recently  Ku  [12],  Dietrich  [5],  and  Deming  [13]  have 
written  more  comprehensive  treatises  on  experimental  error  and  provided  a 
sound  theoretical  basis  for  their  formulations.  Of  these  modern  treatments  Ku 
has  written  the  most  complete  and  thorough  exposition  of  the  propagation  of 
errors.  In  his  paper,  Ku  states  two  forms  the  analysis  of  random  functions 
can  take:  (1)  determining  the  statistical  parameters  of  a given  function  of  a 

random  variable,  or  (2)  determining  the  statistical  parameters  in  which  the 
function  is  assumed  to  tend,  asymptotically  with  a large  number  of 
observations,  to  the  normal  distribution.  Ku's  paper  deals  only  with  the 
second  formulation  of  the  problem.  The  problem,  where  the  function  or  its 
expected  variation  does  not,  for  whatever  reason,  conform  to  asymptotic 
normality  is  of  equal  importance,  however,  and  should  be  addressed.  Section 
5.9  will  address  the  various  methods  which  can  be  used  where  the  assumption  of 
asymptotic  normality  is  not  essential. 

The  basis  for  the  traditional  approach  to  the  study  of  errors  is  embodied 
in  a theorem,  which  is  paraphrased  from  Ku,  as  follows: 

If,  in  some  neighborhood  a function  of  random  variables  is  continuous  and 
possesses  continuous  first  and  second  derivatives  (i.e.  f(x)  eC2),  then  the 
statistical  parameters  of  the  function  can  be  approximated  by 
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f(x)  = expectation  of  f(x)  = f(x) 


(3.8) 


and  the  variance  is  given  by 


1 

var  f(x)  = n 


(3.9) 


(where  x is  a vector  of  n elements)  based  on  asymptotic  normality.  A Taylor 
series  expansion,  from  which  second  and  higher  order  terms  are  dropped,  of  the 
function  over  the  continuous  neighborhood  serves  a basis  for  the  derivation  of 
the  statistical  approximations  to  the  mean  and  variance,  as  given  by  equation 
(3.8)  and  (3.9).  There  are,  then,  three  conditions  attendant  on  the 
acceptance  of  this  development  (1)  that  the  variations  of  the  independent 
variables  are  small  enough  to  justify  ignoring  second  and  higher  order  terms 
from  the  Taylor  series  expansion,  (2)  that  the  approximations  of  statistical 
parameters  tend  to  the  normal  as  the  number  of  observations  becomes  large,  and 
(3)  that  the  function  is  necessarily  well  behaved  in  the  domain  of  interest. 
The  case  of  large  error  or  small  sample  size  is  not  amenable  to  the  above 
conditions,  and  therefore  requires  its  own  methodology.  The  use  of  any  of  the 
familiar  approximations,  which  measures  error  or  computes  a statistical 
parameter,  introduces  an  error  (induced  qualitative  error^)  the  size  of  which 
depends  on  the  formula  employed,  the  number  of  observations,  and  the  magnitude 
of  the  uncertainties  in  the  independent  variables  of  the  function.  The  main 
references,  [ 5 J , [12],  [13],  and  the  subsidiary  references  [14-20]  provide  a 
complete  statement  of  the  theory  of  small  errors.  Of  particular  interest  is 
Ku's  extension  of  the  theory  to  accuracies  of  the  various  stages  of 
approximations . 

^See  Section  2.2.4 
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3. 3. 1.4  Differential  Formulation  of  Error 

An  alternative  formulation  to  the  theory  given,  above  is  to  require  the 
function  to  be  only  of  function  class  Cl , i.e.,  needing  only  a continuous 
first  derivative,  and  expand  the  function  as  a total  differential,  and  use 
this  result  as  a measure  of  absolute  sensitivity.  Let  f(x) , where  x is  a 
vector,  be  , then 

n 

df  = l dXi  (3.10) 

i=l  3xjl 

where  the  x^  are  independent  or  are  functions  of  another  variable.  If  the 
function  is  defined  as 

w = f(u,v)  (3.11) 

where  u is  a vector  of  control  variables  and  v a vector  of  state  variables, 
and  u^  is  one  of  the  u(u^  e U) , the  sensitivity!  to  u^  may  be  developed  as 
follows : 

n m 

dw  = £ du  j + £ dvfc  (3.12) 

j=l  3uj  h=l  3vft 

and  forming  dw/du-^,  we  obtain  (3.13) 

n m 

dw  = a_f  + J duj_  + £ 3_f  dvjj.  (3.14) 

du^  3u^  j = j 3uj  du^  ^=1  3v^  du^ 

as  the  sensitivity  of  the  function  w to  variation  in  u^.  The  first  order 
partial  derivatives  in  equation  (3.14)  are  called  sensitivity  coefficients  in 
the  literature.  The  rates  duj/du^  express  the  functional  dependencies  among 
the  control  elements  on  u^.  The  rates  dv^/du^  express  the  functional 
relationships,  if  any,  of  the  state  variables  with  the  u^  control  variable. 

lln  some  models  it  may  be  important  to  identify  control  and  state  variables 
in  order  to  refine  the  evaluation  of  a sensitivity  analysis.  In  this  example 
u and  v represent  this  distinction. 
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If,  in  equation  (3.14)  we  are  interested  in  time  dependent  changes  in  w,  the 
equation  becomes: 

• 0 m • • 

w - dw  - £ 3 f du-j  + £ 3_f  dvy ; where  “ui  - u^  and  °vk  - v^  . (3.15) 

dt  i=l  3ujL  dt  k=l  3v^  dt  dt  dt 

3. 3. 1.5  Variational  Formulation  of  Errors 

In  analogy  with  the  definition  of  the  differential,  we  define  5 as  the 

variation  in  a function  or  variable.  From  Hildebrand  [21],  the  differential 

of  a function  is  an  approximation  to  the  change  in  that  function  along  a 

particular  curve,  while  the  variation  is  a first  order  approximation  to  the 

d_ 

change  from  curve  to  curve.  The  operators  dx  and  6 are  commutative  for 
independent  variables,  and  are  therefore  interchangeable,  so  that 


— (<$y)  = <5  f * 
dx  \ dx  J 

We  can  write  the  variational  form  of  equation  (3.10)  as 


The  variation  is  usually  defined  as 

fix  = e<J>(x)  (3.22) 

which  changes  a function,  say  y(x) , to  a new  function,  z(x)  by: 

z(x)  = y(x)  + fix  = y(x)  + e<j>(x),  (3.23) 

where  e is  some  appropriate  value. 

Hence  the  values  which  fix  may  take  are  not  restricted.  The  statistical 
approximation  to  the  variance,  as  developed  above  in  Section  3. 3. 1.3,  equation 
(3.10),  can  be  derived  from  (3.21)  by  squaring  6f,  taking  the  expected  values 
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of  both  sides  of  the  equation  and  reducing  the  terms  to  their  variance.  Weld 
[14]  has  derived  the  more  usual  expression  for  the  standard  deviation: 


by  this  method.  For  this  equation  to  be  true,  however,  would  require 
independence  of  the  variables,  as  was  required  in  the  derivation  based  on 
Taylor  series  expansion.  The  Sx  would  be  restricted  to  represent  the  error  or 
variation  in  the  observations,  and  would,  under  the  protocol  of  a large 
sample,  comply  with  the  requirements  of  convergence  to  normality.  This 
condition  is  essential  o all  of  the  formulations  of  classical  error 
analysis . 


2 1/2 


S 


(3.24) 
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4.0  STATISTICAL  METHODS  OF  SENSITIVITY  ANALYSIS 


4.1  Introduction 

This  Section  presents  certain  specialized  applications  and  developments 
of  statistical  methodology  which  are  used  in  the  sensitivity  analysis  and 
evaluation  of  large  models. 

The  methods  discussed  are  sampling  techniques,  which  are  motivated  by 
model  evaluation  problems;  analysis  of  range!  data,  which  is  applicable  to 
parameter  estimation  and  to  the  evaluation  of  large  errors;  the  use  of  partial 
rank  correlation  coefficients  as  measures  of  sensitivity;  and  three  standard 
statistical  methods,  range!  analysis,  minimum  variance,  and  regression,  which 
are  used  for  estimating  and  sizing  sensitivity. 

Of  these  subjects  the  sampling  techniques  and  application  of  partial  rank 
correlation  coefficients  as  measures  of  sensitivity  are  of  particular  value  to 
those  interested  in  assessing  a large,  complex  computerized  model.  The  best 
of  the  sampling  schemes,  latin  hypercube  sampling,  is  a procedure  which 
structures  input  trials  in  an  optimum  and  comprehensive  model-testing  design. 

4.2  General  Sampling  Techniques 

Three  principal  methods  of  sampling  are  presented.  Two  of  the  three 
methods  are  familiar  and  are  the  subjects  of  a wide  and  comprehensive 
literature.  The  third  method,  however,  is  comparatively  new,  that  of  latin 
hypercube  sampling  (LHS) , and  is  discussed  and  compared  with  standard  sampling 
procedures . 


1- 


range"  is  used  here  to  refer  to  the  extreme  values  of  an  interval  around 
value. 


a 
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4.2.1  Random  Sampling 


In  random  sampling,  values  for  input  variables  are  determined  by  a scheme 
which  assigns  a value  equally  likely  over  all  the  values  the  input  variable  is 
permitted  to  have.  Further,  the  sample  values  so  selected  have  the  same 
probability  of  selection  as  any  other  set  of  assignments  of  values  to  the 
variables.  The  collection  of  values  over  all  the  variables  is  called  the 
sample  space.  The  sample  then  is  the  set  of  assigned  values  to  each  input 
variable  under  the  condition  that  each  assignment  was  as  equally  likely  as  any 
other  assignment  and  that  each  assignment  is  independent  from  all  other 
assignments . 

4. 2. 1.1  Evaluation  There  exists  an  immense  literature  on  the  analysis  of 
random  sampling  and  its  use  in  estimating  the  statistical  parameters  of  the 
population  from  which  the  sample  was  taken.  Our  intent  in  introducing  random 
sampling  is  to  compare  its  efficiency  against  other  methods  of  sampling, 
particularly  as  it  relates  to  generating  data  for  senstivity  analysis  and 
model  evaluation.  It  will  be  shown  below  that  a variance  based  on  a random 
sample  will  be  greater  than  or  equal  to  the  variance  based  on  a stratified 
sampling  scheme.  See  references  [22]  and  [23]  for  technical  exposition. 

4. 2. 1.2  Application  The  use  of  random  sampling  is  wide  spread,  and  in  the 
specified  requirements  of  some  models  it  has  the  desired  properties  needed  to 
attain  certain  conditions,  but  as  a technique  for  generating  data  for  model 
evaluation  it  is  the  most  costly  and  inefficient  of  the  three  methods  outlined 
in  this  Section. 
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4.2. 2 St  ra  t i ft  ed  Samp  1 i ng 

In  stratified  sampling  < lie  sample  space  of  input  variables  is  partitioned 

into  n subspaces.  Following  (lie  development  by  McKay,  Beckman,  and  Conover 

I 23),  the  sample  space  is  partitioned  into  disjoint  sets  Sj  of  size  Pj , and 

if  Xj  j are  a random  sample  (j)  on  Sj , and  with  y = h(x),  an  unknown  but 

observable  mapping,  then  the  strata  means  and  variances  of  g(y)  are: 

Pi  = F(tf(yn))  = / g(y)  }_  f(x)dx  (4.1) 

Si  pt 

«t2  = V(  g(  y i i ) ) = / (g(y)“Pi>2  f ( x)dx  (4.2) 

Si  Pi 

and  the  estimator  over  the  sample  space 

N 

T = J_  l g(Ul)  , (4.3) 

N i=  1 


where  g (uj.)  is  an  arbitrary,  known  function.  Then 


(4.4) 


and  the  variance  of  Ts  is  given  by: 

V(TS)  = l / £±1  1 °i2  , (4.5) 

i=l  | ni  | 

where  the  "s"  on  T as  a subscript  indicates  "stratified  sampling"  and  I is  the 
number  of  disjoint  subsets  S^.  Invoking  the  development  of  Tocher  [24], 
equation  (4.5)  becomes: 

I 

V(TS)  = V(Tr)  - £ I Pi  (Pi-T)2  (4.6) 

N i=l 

where  "R"  on  TR  denotes  the  estimator  based  on  random  sampling  and  x is  the 
estimate  of  the  mean  under  stratified  sampling.  It  follows  from  equation 
(4.6)  that 
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< 

V(TS)  = V(Tr) 


(4.7) 


which  is  the  desired  demonstration  that  stratified  sampling  is  a more 
efficient  method  than  is  random  sampling.  A necessary  condition  attendant  to 
this  conclusion  is  that  the  sample  sizes,  n^ , are  chosen  so  that 

n±  = P-jN  (4.8) 

which  provides  for  proportional  allocation.  The  variance  reduction  shown  in 
equation  (4.6)  is  a function  of  the  differences  between  the  strata  means  (y^) 
and  the  overall  mean  (t). 

4.2.2. 1 Evaluation  Although  the  advantage  of  a reduced  variance  is  gained 
through  the  use  of  stratified  sampling  it  is  not  certain  that  this  method  of 
sampling  can  be  applied  to  large  simulations.  The  sample  size  per  cell, 
equation  (4.8),  and  the  partitioning  of  the  data  set  may  be  impossible  to 
impose  or  may  not  be  a natural  realization  of  the  system  or  subject  which  has 
been  modeled.  The  amount  of  the  reduction  is  also  a factor  to  consider.  In 
both  the  random  and  stratified  sampling  methods  N runs  are  made.  If  the  cost 
of  partitioning  is  substantial  it  must  be  measured  against  the  value  attained 
in  the  reduction  of  the  variance.  It  follows,  then,  that  a cost-benefit 
analysis  imposes  on  the  decision  to  use  stratified  sampling. 

4.2.3  Three  Point  Sampling 

This  plan  uses  three  values  to  represent  the  range  of  values  for  inputs. 
These  values  are  assigned  independently  to  each  run  based  on  a probability  for 
each  value,  the  relationship  being 

P1  = p2  = I (l-po>  (4.9) 

2 
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where  P0  is  the  probability  of  the  median  value  which  is  assigned  by  the 

analyst.  Under  proportionate  allocations  3(l-p0)_^  values  are  assigned  at 

2 

random  to  n runs  of  the  model.  For  example  if  x_ , x+,  and  xQ  are  the  lower, 
higher  and  median  values  of  variable  x with  probabilities  Pj , P2 , and  P0 , then 
a random  draw  based  on  an  interval  in  proportion  to  the  values  of  , P2 , and 
PQ  will  determine  which  of  the  three  values  will  be  selected  for  a given  run. 

This  procedure  has  been  used  to  evaluate  the  COAL  2 National  Energy 
Model,  and  to  the  author's  knowledge,  has  not  been  used  except  in  this  one 
instance.  The  source  for  this  description  is  Ford,  Moore,  and  McKay  [25]. 
4.2.3. 1 Evaluation  The  three-point  sampling  plan  is  actually  only  a special 
case  of  stochastic,  factorial,  stratif ied-sampling  methodology.  In  this 
technique  each  input  variable  is  assigned  a density  function,  say  f^(x) , for 
x-^L  < x < x^H.  Values  of  x in  this  range  are  defined  such  that  the  following 
relation  is  obtained: 

xiL  = aio  < ail  < < aik  = XiH  (4.10) 

and  the  probability  of  the  interval  (a^  j_j , ai,j)  is  1/k.  These 
intervals,  one  for  each  input  variable,  are  sampled  to  produce  different 
values  of  input  variables  for  each  computer  run. 

There  are  two  major  objections  to  this  sampling  method.  The  first  is  the 
practical  difficulty  of  assigning  a density  function  for  each  input  variable, 
particularly  for  a model  with  a large  input  data  base.  The  second  objection 
arises  in  the  possible  re-selection  or  assignment  of  values  which  may  have 
been  used  in  a prior  run.  This  is  the  problem  of  sampling  with  replacement. 
This  problem,  along  with  the  realization  that  the  sampling  procedure  may  not , 
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based  on  probability  selection,  provide  a proper  spread  of  values  which  are 
needed  to  evaluate  the  model  performance,  creates  the  possibility  that  some 
degree  of  monitoring  of  values  assigned  is  necessary  and  also  that  the  number 
of  needed  runs  may  actually  be  increased  beyond  initial  expectations. 

The  next  Section  will  discuss  how  these  difficulties  are  precluded  in  the 
Latin  Hypercube  Sampling  method. 

4. 2. 3. 2 Application  This  scheme  of  sampling  is  for  statistically  independent 
input  variables  and  models  for  which  the  preparation  of  the  scheme  is  not 
excessively  costly  in  time  and  resources.  The  number  of  runs  for  this  type  of 
statistical  analysis  is  of  the  form 

N * k1  (4.11) 

where  k is  the  number  of  intervals  and  I the  number  of  input  variables.  For 
some  models,  those  of  long  running  times,  a large  N may  prohibit  the  use  of 
this  approach. 

4.2.4  Latin  Hypercube  Sampling  Method 

The  principal  references  on  this  method  are  McKay,  Conover,  and  Whitman 
[22],  and  McKay,  Beckman,  and  Conover  [23]. 

When  using  the  stratifed  sampling  methods  there  are  many  options  on  the 
choice  of  the  number  of  intervals  (k)  of  each  input  variable  and  the  fraction 
of  an  experimental  design  wherein  only  some  of  the  combinations  of  interest 
are  actually  evaluated.  The  Latin  Hypercube  Sampling  (LHS)  method  uses  the 
same  number  of  intervals  for  each  input  variable  as  the  number  of  computer 
runs  to  be  made.  The  intervals  for  each  input  are  assigned  by  defining  a 
sequence  of  independent,  uniform  random  variables,  ranking  them  from  lowest  to 
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highest,  and  then  use  this  ranking  to  determine  the  computer  run  in 

which  that  input  value  will  be  used.  As  in  the  stratified  sampling  method  the 

intervals  must  be  sampled  to  obtain  specific  values  of  the  input  variables. 

The  values,  however,  may  be  fixed  by  the  model  user  instead  of  having  the 

assignment  dependent  on  an  interval  density  function.  This  choice  is  an 

important  option  in  terms  of  cost  and  preparation.  This  method  can  be 

regarded  as  an  extension  of  Latin  Square  sampling. 

Reference  [23]  demonstrates  that  the  variance  of  an  estimator  based  on 

LHS,  in  terms  of  the  variance  based  on  random  sampling,  is: 

V(Tl)  = V(Tr)  + (N-l)  (N_k(N-l)"k)  l (Pi-T)  (prT)  , (4.12) 

N 

where  £ is  over  the  restricted  space  of  pairs  (p^.pj)  having  no  coordinates 

< 

in  common,  and  V(Tl)  = V(Tr)  if  and  only  if 

N-k(N-irkJ  (ui-T)  (p-j-T)  = 0 (4.13) 

No  direct  way  of  comparing  the  variance  of  an  estimator  based  on  the  LHS  to 
that  based  on  stratified  sampling  is  available.  However  Reference  [23] 
states,  and  proves  the  assertion,  that  if  a function  Y = h(x^,  ••*,  xR)  is 
monotonic  in  each  of  its  arguments,  and  g (Y)  (see  equation  (4.4))  is  a 
monotonic  function  of  Y,  then 

< 

V(Tl)  = V(Tr)  . (4.14) 

4. 2. 4.1  Evaluation  The  principal  advantage  of  LHS  is  to  overcome  the 
difficulties  of  the  stratified  sampling  method  by  insuring  a representative 
coverage  of  all  input  values  within  the  set  of  prescribed  runs.  The  output  of 
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the  sampling  design  uses  all  of  the  intervals  on  the  range  of  each  input 
variable.  The  proofs  and  mathematical  developments  given  in  References  [22] 
and  [23]  assume  a random  selection,  without  replacement,  for  an  input  value. 

It  appears,  however,  that  a predetermined  assignment  of  values  from  a 
given  partition  on  each  input  interval  could  be  made  which  would  allow  a 
direct  study  of  particular  effects  associated  with  certain  values  of  the  input 
data.  This  technique  uses  LHS  as  a probe  for  special-case  examinations,  by 
avoiding  the  necessity  to  extract  such  information  from  a more  general 
statistical  design.  This  application  fixes  the  number  of  runs,  however,  and  a 
tradeoff  between  the  number  of  runs  (or  model  evaluations)  and  the  interval 
partition  must  be  evaluated.  To  assign  values  might  violate  statistical 
requirements,  but  it  would  provide  limited  but  important  information  on  model 
performance . 

The  advantage  of  the  LHS  is  revealed  when  the  model  output  is  dominated 
by  only  a few  input  variables.  But  the  LHS  ensures  that  each  component  is 
fully  stratified,  no  matter  which  variables  may  or  may  not  dominate  the 
results . 

In  Reference  [22]  McKay,  et  al , conclude  that  factorial  stratified 
sampling  (FSS)  and  Latin  hypercube  sampling  are  clearly  superior  to  estimates 
based  on  random  sampling.  The  amount  of  gain,  however,  depends  on  the  model 
under  assessment,  but  in  their  experience  the  LHS  resulted  in  a maximum 
variance  reduction  of  about  50  to  1 over  random  sampling.  Early  studies  also 
indicate  LHS  is  "far  superior"  to  FSS. 

4.3  Analysis  of  Range 

In  a basic  paper  Tsukibayashi  [26]  has  developed  the  moments  of  the  range 
for  rectangular,  triangular,  exponential,  gamma,  and  normal  distributions 
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His  results  include  a calculation  of  certain  coefficients  to  be  used  in 
estimating  the  first  four  moments  of  the  distribution. 

Tsukibayashi  states  that  the  moments  of  the  range  are  functions  only  of 
the  sample  size  N and  the  variance  a;  the  following  are  derived.  If  Wjj  is  the 
range,  then: 

Pj  = E(Wn)  = djgO  , and  (4.15) 

1*2  = E(Wjj2)  = Cjj2a2  ; (4.16) 

higher  moments  are  not  given  here.  The  table  below  summarizes  the  results  of 
the  analysis  of  range  in  terms  of  statistical  parameters. 


Distribution 

<*N 

CNZ 

Rectangular,  f(X)=l ,0<X<1 ,o^= 1/12 

12  ( n-1\ 

12Pln-^  1 

\ trfl/ 

|^(n+l)(n+2)  J 

Triangular  , f (X)=2X  ,0<X<  1 ,a2=l  / 1 8 

1*  |~2n  - lB(l/2,n+l)| 

1 

t 

CM 

1 

1 

00 
*— H 

(B  is  the  beta  function) 

2n+l  2 J 

2(  rt+1 ) 1 

B ( 1 / 2 , rH-1 ) I 

n-1 

n-1 

Exponential , f (X)  = e-X  ,0<x<°°  ,a2=l 

I 1 

V+I  i, 

k=l  k 

k=l  k2 

For  the  more  complicated  expressions  for  the  gamma  and  normal  distributions  of 
the  range  the  reader  should  refer  to  [26],  pages  66  and  67.  The  theory* 
underlying  the  coefficients  given  above  was  developed  by  Tippett  and  Pearson. 
4.4  Minimum  Variance  Sampling 

For  models  which  simulate  the  states  of  a system  or  activity,  such  as 
Markov  processes,  queues,  or  stochastic  processes,  the  sampling  is  usually 
performed  in  proportion  to  their  natural  freauency  of  occurrence.  A technique 
has  been  developed  by  Bayes  [27]  which  allows  the  sampling  frequency  of  the 

* The  reference  to  Tippett  is  Biometrika,  17  (1925). 
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states  to  be  independent  of  their  natural  frequency.  A sampling  theory  of 
Markov  chains  is  given  by  Bayes  which  allows  some  statistics  of  the  state 
frequencies  to  be  estimated  with  minimum  variance  for  a given  sample  size. 

In  models  of  this  kind  some  states  are  experienced  infrequently,  which 
requires  a large  number  of  runs  to  collect  enough  data  on  these  states  to 
perform  a proper  evaluation.  This  situation  has  led  to  the  development  of 
"importance  sampling"  in  which  the  simulation  is  constrained  to  experience 
rare  states  independently  of  the  expected  frequency  of  the  number  of  times 
these  states  are  actually  realized  in  a normal  use  of  the  simulation. 

Bayes  applies  the  theory  by  approximating  the  model,  in  this  case  a 
queueing  problem  within  a teleprocess  simulation,  by  a Markov  chain  in  which 
states  are  assigned,  transition  probabilities  are  given,  and  means  and 
variances  of  intervals  spent  in  each  state  are  assigned.  The  statistical 
properties  of  the  Markov  chain  surrogate  model  so  defined  will  be  similar  to 
those  of  the  original  simulation. 

4.5  Regression 

The  subjects  of  linear  and  nonlinear  regression,  and  polynomial  curve 
fitting  are  so  familiar  that  the  structures  of  these  techniques  will  not  be 
repeated  here.  This  Section  is  included  for  the  purpose  of  reviewing  the 
contributions  of  regression  to  the  sensitivity  analysis  of  models. 

A proper  regression  is  known  to  reduce  the  length  of  a confidence 
interval  and  to  influence  the  value  of  the  model  output,  assuming  a causal 
relationship.  The  regression  model  itself  contributes  information  on 
sensitivity  through  the  magnitudes  of  the  coefficients  and  the  numbers  of 
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variables  taken  into  the  regression  formulation.  Both  curve-fitting  and 
regression  provide  functional  relationships  which  can  be  used  according  to  the 
theory,  given  above,  on  error  propagation  and  analysis. 

We  quote  an  important  comment  from  Deming  [13]  on  the  matter  of  indirect 
contributions  to  knowledge  based  on  regression  analysis:  "There  have  been 

many  instances  when  deductions  made  from  a fitted  curve,  or  from  a series  of 
curves,  have  made  it  unnecessary  to  perform  certain  other  experiments.  As  an 
instance  ...  where  a quartic  was  fitted  to  compressibility  data  on  carbon 
dioxide  ...  this  quartic  gave  data  on  the  index  of  refraction,  the 
Joule-Thomson  coefficient,  entropy,  and  other  physical  properties  that  would 
be  difficult  and  time  consuming  if  direct  observations  were  required." 

4.6  Partial  Rank  Correlation  Coefficient  (PRCC) 

McKay,  et  al  [22],  have  proposed  and  used  the  PRCC  as  an  alternative  to 
the  derivative  approach.  The  partial  correlation  coefficient,  originally 
employed  to  measure  linear  relationships  among  variables,  is  adjusted  to  use 
ranked  data,  and  allows  for  nonparametric  tests.  Following  [22]  the 
formulation  is  as  follows: 

Let  (xj^,  X2i,  ...,  x^,  y^)  denote  the  values  of  k input  variables 
and  the  output  variable  for  the  ith  computer  run,  i=l,  ...,  n.  Let 
(xrli,  xr2i,  ...,  xrki,  yr^)  denote  the  associated  rank  variables. 

For  example,  if  n=3,  and  xn=12,  x12=25,  x13=10,  then  xr11=2,  xr12=3,  and 
xri3=l. 

Let  C be  a (k+l)x(k+l)  symmetric  matrix  of  elements  cjj  , when 
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I,  <xrit-*ri><xrjt-*rj> 
:ij  ' 


n _ n _ “1  1/2 

I (Xrit-*ri>2  l (X^s-x/)2 

t=l  S=1 


for  i,  j = 1 , 2,  ...»  k, 


and 


(4.17) 


ci >k+l 


JE  (xrn-*ri)(yrt-yr) 


I ttrit-*ri>2  l <yrs-yr>2 

t=l  S=1 


1/2 


(4.18) 


where  c^^=l  and  Cj:j=Cji  for  i ,j=l , . . . ,k+l . 

The  values  xr^  and  yr  are: 

_ n _ n 

xri=£  l = n+1  , yr  = £ l yrt  - fli  • <4-19> 

n t=  1 2 n t=  1 2 

Let  B=[b^j]=C_l  then  the  PRCC  between  input  variable  i and  output  variable  y 
is  given  by: 

rPiy  = -bi>k+1  ( bi;L  bk+1  ,k+1)-1/2  (4.20) 

The  principal  reference  given  in  [22]  is  Steel  [29]  for  the  underlying 
theory  for  this  approach.  References  [28]  and  [30]  provide  detailed  accounts 
of  the  statistical  theory  of  partial  correlation  coefficients.  Reference  [30] 
describes  the  theory  from  the  perspective  of  geometric  constructions  in  which 
all  the  relations  are  expressed  in  terms  of  lengths  of  vectors  and 
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subtended  angles,  thereby  portraying  the  theory  of  partial  correlation  with 


the  trigonometry  of  an  n-dimensional  constellation  of  points. 
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5.0  ANALYTICAL  METHODS  OF  SENSITIVITY  ANALYSIS 


5.1  Introduction 

This  section  describes  methods  of  sensitivity  analysis  which  are 
analytical  formulations  rather  than  statistical  designs,  and  are  based  on 
special  needs  that  traditional  methodology  could  not  meet.  In  some  cases  the 
methods  are  experimental  and  of  a singular  nature,  in  others  of  a more  solid 
nature  which  advance  the  state  of  the  art  through  the  introduction  of 
innovative  techniques  or  by  the  further  development  of  established  sensitivity 
analysis. 

5.2  Criterion  Sensitivity  Analysis 
5.2.1  Descripton 

This  method  has  developed  at  the  MIT  Energy  Laboratory,  Model  Assessment 
Group.  The  single  source  is  Schwippe  and  Grubl  [31],  and  seems  not  to  have 
had  much  of  a life  after  its  genesis  in  the  mid-^TO's.  In  [31]  the  authors 
describe  the  methods  as  follow: 

The  model  is  be  represented  as  y=f (x) , where  x is  a vector  of 
exogenous  input  parameters  and  y is  a vector  of  model  outputs.  The  perturbed 
output  is  y(Ax)  = f(x+Ax) , where  Ax  is  the  perturbation,  y(Ax)  is  the 
perturbed  output,  and  6y(Ax)  = y(Ax)-y(x)  is  the  output  perturbation.  Assume 
that  some  particular  scalar  criterion  function  C of  the  output 
perturbation 

C(Ax)  = C[ <5y( Ax)  ] (5.1) 

has  been  selected.  Then  the  basic  idea  of  the  criterion  sensitivity  analysis 
is  to  find  that  Ax  which  maximizes  C(Ax).  This  can  be  viewed  as  a "worst 
case"  analysis  which  finds  the  maximum  sensitivity  as  measured  by  C(Ax).  The 
MIT  approach  was  to  choose  C(Ax)  from  one  of  the  perturbed  outputs  6y(Ax). 
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The  authors  emphasize  the  usefulness  of  this  method  for  models  of  a highly 
nonlinear  nature. 

5.2.2  Evaluation 

This  method  was  developed  because  simple  linearization  is  not  appropriate 
for  many  models.  Criterion  sensitivity  analysis  addresses  in  particular  the 
problem  of  nonlinearity  in  the  function-model  represented  by  equation  (5.1). 

An  objection  noticed  by  the  authors  is  the  difficulty  in  characterizing  the 
input  uncertainties  (i.e.  defining  the  space  of  AX  or  its  probability 
density).  The  advantage  of  the  method  is  that  it  provides  explicit  mechanics 
for  dealing  with  input  uncertainties  in  a systematic  fashion.  No  further 
descriptions  of  this  method  were  available  beyond  [31]  at  the  time  this  survey 
was  prepared. 

5.3  Describing  Function  Sensitivity  Analysis 
5.3.1  Description 

Our  source  for  the  description  of  this  technique  is  the  same  as  for  the 
criterion  method,  Reference  [31],  which  is  the  basis  for  the  following 
description: 

With  the  describing  function  approach,  the  input  perturbation  vector,  Ax, 
is  taken  as  a random  vector  characterized  by  its  probabilitv  density  p(Ax). 

Any  density  function  is  permitted  and  it  is  not  assumed  that  the  members  of  Ax 
are  independent. 

The  describing  function  is  defined  to  be  a vector  polynomial  function  of 
the  vector  Ax,  as  follows 

M 

D(  Ax)  = Aq+AjAx*  J e£nAxtA2mAx  + ...  (3.2) 

m=l 
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where  M is  the  dimension  and  em  is  a unit  column  vector  (all  zeros  except  1 in 
the  mth  row).  Equation  (5.2)  is  simplifed  to 

D(Ax)  = A<J>(Ax)  (5.3) 

where  <J>(Ax)  is  a vector  of  "l"s,  Ax,  and  powers  and  cross  products  of  Ax  up  to 
the  number  of  terms  desired,  of  dimension  K,  and  A is  a matrix  of  Aq,  Aj , 
etc. 

The  describing  function  problem  is  then  defined  as  follows:  Find  values 

of  matrix  A and  the  number  of  K terms  such  that 

6y(Ax)  « D(Ax)  (5.4) 

Without  reproducing  the  remaining  mathematical  development  the  following 
concluding  paragraph  is  cited  [31,  page  506]  to  indicate  the  detail  and  intent 
of  this  approach. 

"In  summary,  the  preceding  process  involves  fitting  a set  of  describing 
functions  to  a set  of  model  runs  that  represent  a set  of  points  on  the  model's 
input-output  response  surface.  The  suggested  manner  of  determining  the  best 
fit  has  been  a weighted  least  squares  approach,  where  the  importances  of  the 
different  points  are  weighted  by  the  likelihood  that  the  response  is  going  to 
be  in  the  neighborhood  of  those  points.  With  a large  enough  set  of  describing 
functions,  that  is  as  K approaches  N* , the  fitting  of  the  different  points  can 
be  "perfect.”  Such  "perfect"  fits  are  highly  susceptible  to  capitalization  on 
chance  effects,  and  the  preceding  discussions  suggests  that  one  way  of 
avoiding  such  spurious  fits  is  to  restrict  the  number  of  describing  functions 

In  is  the  number  of  perturbations,  taken  larger  than  K. 
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in  the  set  being  used.  There  are  endless  variations  on  this  particular 
suggested  procedure  for  determining  the  "best"  fit.  A few  of  these  variations 
include : 

(1)  minimizing  the  maximum  residuals  between  the  surface  and  the  points, 

(2)  minimizing  absolute  differences  or  some  other  robust  measure  rather 
than  least  squares,  or 

(3)  using  weighted  maximum  residuals  or  weighted  robust  techniques. 

In  addition  it  would  be  possible  to  exploit  any  intuitions  one  might  have 
about  the  true  characteristics  of  the  model's  response  surface.  For  example, 
if  a large  number  of  points  were  available  and  one's  intuition  suggested  that 
the  response  surface  should  be  a relatively  smooth  connection  of  those  points, 
then  the  fitting  criterion  might  be  to  minimize  the  integral  of  the  deviation 
between  the  fitted  surface  and  the  piecewise  linear  connection  of  the 
available  points  with  their  nearest  neighbors  (or  the  supporting  hyperplane 
n-gons  connected  over  all  convex  sets  of  available  points).  The  principal 
drawback  to  these  more  elegant  techniques  is  that  they  may  not  be  nearly  as 
easily  solved  as  is  the  least  squares  approach." 

5.3.2  Evaluation 

In  addition  to  difficulties  of  establishing  the  character  of  the  input 
uncertainties  and  the  final  comment  of  the  cited  paragraph  given  above,  the 
authors  assert  that  describing  functions  can  be  very  useful  in  cases  where  the 
direction  and  magnitude  of  the  changes  in  outputs  is  of  greater  interest  than 
the  relative  effects  of  different  types  of  hypothesized  inputs. 
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5.4  Adjoint  Sensitivity  Analysis 
5.4.1  Description 


A recent  and  important  method  of  sensitivity  analysis  has  emerged  from 
the  collective  work  of  scientists  and  model  builders  at  the  Oak  Ridge  National 
Laboratory  (ORNL).  This  method  is  called  adjoint  sensitivity  analysis  and  had 
its  origin  in  some  earlier  papers  by  Oblow  [33,34]  and  Cacuci  [36]  which  dealt 
with  general  sensitivity  analysis  of  problems  in  nuclear  cross-section, 
shielding,  and  reactor  physics. 

The  theoretical  method  was  recognized  for  its  potential  as  a practical 
tool  for  evaluating  the  large  economic  - energy  models  which  were  also  in  use 
at  ORNL. 

The  basic  paper  from  which  our  description  is  taken  is  Alsmiller,  et  al , 
[32].  It  provides  a clear  statement  of  the  essentials  of  the  method.  The 
underlying  theory  and  mathematics  are  given  in  [33-37],  Reference  [38] 
describes  the  application  of  the  method  to  an  ORNL  model  called  LEAP. 

The  following  text,  which  describes  adjoint  sensitivity  theory,  is  taken, 
with  minor  editing  for  clarity,  from  Alsmiller  [32],  pages  3-7. 

Let  the  system  under  consideration  be  represented  by  N nonlinear 
algebraic  equations  in  N unknowns  which  may  be  written  symbolically  as 

F(p,x)  = 0,  (II. 1) 

where  F is  an  N-vector  function  of  the  unknown  dependent  variable  N-vector  p 
and  all  data  elements  are  represented  by  the  vector  x.  It  is  to  be  understood 
that  the  number  of  elements  in  x is  not  at  all  related  to  N and  can  be 
substantially  larger.  In  the  following  it  is  assumed  that  for  a specific 
choice  of  x a unique  solution  of  Eq.  (II. 1)  exists  and  is  represented  by  p. 

The  vector  p is  thus  a function  of  x. 
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Let  R represent  any  result,  hereinafter  also  called  a response,  of  the 

calculations  that  are  of  interest.  In  particular,  let 

R = R(p,x),  (II. 2) 

wehre  R(p,x)  is  any  given  nonlinear  function  of  p and  x.  The  response  R is  a 

scalar  that  may  be  calculated  from  Eq.  (II. 2)  when  the  solution  p of  Eq . 

( II . 1 ) has  been  obtained  for  a given  specification  of  x. 

If  now  x is  defined  to  be  a general  element  of  x the  sensitivity  of  R to 

x is  defined  to  be  and  the  general  problem  of  sensitivity  theory  is  to 

dx 

determine  for  each  and  every  one  of  the  elements  of  x.  To  accomplish 
dx 

this,  differentiate  Eq.  (II. 2)  with  respect  to  x,  to  obtain: 

dR=  9E+9R#  dp  „ Q . (II. 3) 

dx  9x  9p  dx 

where  is  a row  vector.  The  first  term  on  the  right-hand  side  of  of  Fq. 

9p 

(II. 3)  is  called  the  "direct  effect"  of  x on  R and  the  second  term  on  the 

right-hand  side  in  Eq.  (II. 3)  is  called  the  "indirect  effect"  of  x on  R,  since 

it  reflects  the  implicit  effect  of  x on  R through  p. 

The  evaluation  of  the  second  term  on  the  right-hand  side  of  Eq.  (II. 3) 

begins  with  the  differentiation  of  Eq.  (II. 1)  with  respect  to  x to  give 

Hi  + = 0.  (II. 4) 

9x  9p  dx 

Now  defining  two  new  N-vectors 
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(II. 5) 


$ = dp 
dx 

s H _ 9F  , (II. 6) 

3x 

and  the  N by  N matrix 

A E 3F  , (II. 7) 

3p 

Eq.  (II. 4)  may  be  rewritten  as 

A$  = S . (II. 8) 

Here  $ is  the  differential  change  in  the  solution  p with  respect  to  a change 

in  x and  it  solves  the  linear  inhomogeneous  equation  whose  source  is  the 

negative  of  the  partial  differential  change  in  F with  respect  to  x.  The  fact 

that  A is  a linear  operator  is  clear  from  Eq.  (II. 4)  since  it  cannot,  in 

principle,  depend  on  ^_P  . The  matrix  A does  depend  on  both  p and  x. 

dx 

The  expression  for  in  Eq.  (II. 3)  can  now  be  evaluated  by  solving  Eq. 

dx 

(II. 8)  for  $.  The  difficuly  here  lies  in  the  fact  that  S depends  explicitly 

on  the  particular  dx  being  considered  and  thus,  in  general,  very  large  systems 

of  linear  equations  represented  by  Eq . (II. 8)  must  be  solved  for  each  dx  to 

obtain  dR  . For  very  large  systems  where  x contains  thousands  of  elements 
dx 

this  is  impractical.  It  should  be  noted  that  the  explicit  assumption  is  being 
made  here  that  it  is  not  feasible  to  numerically  construct  the  inverse  of  the 
matrix  A.  If  this  inverse  could  be  obtained  then  the  need  for  solving  the 
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large  systems  of  linear  equations  many  times  would  be  ameliorated  and  the 
above  objection  would  not  be  valid. 

The  adjoint  formulation  of  sensitivity  theory  is  of  interest  because  it 
avoids  the  necessity  of  solving  a large  set  of  linear  equations  for  every  x in 
x.  The  set  of  equations  adjoint  to  the  Eq . (II. 8)  will  be  written 


where  a superscript  * is  used  to  denote  adjoint , a "t"  over  a symbol  is 


inner  products  of  with  Eq.  (II. 9)  and  fl**1  with  Eq.  (II. 8)  are  made  one 
obtains,  by  using  Eq . (II. 10): 

*t  • S*t  = • S , 


(II. 9) 


where  the  N by  N matrix  A*  is  defined^  to  be: 


A*  = At, 


(II. 10) 


used  to  indicate  a transpose , and  where  S*  is  yet  to  be  defined.  If  now  the 


and  if  S*  is  defined^  by: 


3RC 

S*  = 3 p 


(11.12) 


then  using  Eqs.  (II. 11)  and  (II. 12),  Eq . (II. 3)  may  be  written 
dR  = 3 R + <j>*t  • g, 


(11.13) 


dx  3x 


where  is  the  solution  of  Eq . (II. 9),  with  S*  given  by  Eq.  (11.12). 


^The  proof  of  this  definition  is  given  in  Cacuci , et  al , [36].  The  matrix  A* 
is  formed  by  transposing  the  adjoints  of  the  components  of  A. 

^The  justification  for  this  definition  is  given  in  Cacuci,  et  al,  [36]. 
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Equation  (II. 13)  is  the  important  result  of  adjoint  sensitivity  theory. 1 

The  quantity  $ must  be  obtained  by  solving  the,  often  large,  set  of  linear 

equations  given  in  Eq.  (II. 9),  but  this  may  be  used  in  Eq.  (11.13) 

for  x of  the  vector  x.  Thus,  one  calculation  of  $ enables  the  calculation 

of  for  all  x in  x for  a specific  R.  Each  new  response,  R,  that  is 
dx 

considered  requires  a new  calculation  because  of  Eq.  (11.12). 

5. A. 2 Evaluation 

The  observations  on  the  effectiveness  and  efficiency  of  adjoint 
sensitivity  analysis  are  taken  from  Alsmiller  [38]  , with  some  minor 
emendations  and  omissions. 

Several  important  observations  should  be  made  in  regard  to  the  results. 
The  first  advantage  of  this  approach  is  that  the  adjoint  equation  is 
independent  of  any  operation  involving  differentiation  with  respect  to  x. 

This  property  means  that,  no  matter  how  large  the  set  of  input 

parameters  is,  only  a single  adjoint  equation  needs  to  be  solved  to  compute 

dR/dx  for  all  x in  x. 

This  is  in  contrast  to  any  direct  method  of  computing  dR/dx  by  changing 
each  x by  a finite  amount  (i.e..  Ax)  and  solving  Eq.  (II. 1)  to  approximately 
get  dR/dx  = Ar/Ax  or  by  using  Eq.  (II. 8)  directly  to  solve  exactly  for  dp/dx. 
In  both  of  these  latter  cases,  the  solution  of  a large  system  of  equations 
[i.e.,  either  Eq.  (II. 1)  or  II. 8)]  is  needed  for  each  parameter  x to  be 
studied.  The  adjoint  approach  is  therefore  extremely  economical  to  use,  and 

^An  alternative  derivation  of  equation  11.13  is  given  in  Alsmiller,  Jr., 
et  al.  , [53],  This  derivation  is  based  on  an  algebraic  representation  rather 
than  on  linear  operator  theory. 
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with  a very  large  set  of  input  parameters,  it  may  be  the  only  approach  to 
obtaining  all  sensitivity  coefficients. 

The  second  advantage  is  that  the  adjoint  equation  is  a linear  equation  in 
<}>  , which  is  generally  much  easier  to  solve  than  the  repeated  calculations 
using  Eq . (II. 1),  which  is  nonlinear.  The  third  advantage  is  that  this 
approach  allows  the  exact  sensitivity  coefficient  dR/dx  to  be  evaluated  by 
using  Eq.  (11.13)  (i.e.,  no  perturbations  in  x and  no  approximations  such  as 
AR/Ax  = dR/dx  are  needed). 

On  the  other  hand,  certain  limitations  and  disadvantages  of  this  method 
are  also  apparent.  The  first  and  foremost  disadvantage  is  that,  to  construct 
the  adjoint  equations  [i.e.,  Eq . (II. 9)]  and  evaluate  dR/dx  with  Eq.  (11.13), 
various  derivatives  of  f and  R must  be  evaluated  analytically.  This 
evaluation  requires  an  in-depth  knowledge  of  the  functional  form  of  the 
complete  set  of  equations  described  by  f.  Second,  the  structure  of  the  A* 
matrix,  although  independent  of  <j> * , is  not  independent  of  p.  This  means  that 
the  solution  of  Eq . (II. 1)  must  be  available  to  construct  A*,  thus  tying  the 
adjoint  equations  to  the  solution  of  the  forward  equations.  Therefore,  to 
evaluate  Eq . (11.13)  for  dR/dx,  both  solution  vectors  p and  <{>*  must  be 
available  in  addition  to  all  analytic  derivatives  of  f and  R with  respect  to 
both  p and  x in  x.  A tradeoff  study  between  analytic  work  in  the  adjoint 
approach  and  computational  time  in  the  alternative  forward  approaches  may  be 
necessary;  the  usefulness  of  either  method  will  depend  on  the  circumstances  of 
the  sensitivity  study  being  contemplated. 
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One  additional  disadvantage  of  the  adjoint  approach  is  that  an  adjoint 
solution  is  required  for  each  new  response  whose  sensitivity  is  being  studied. 
This  is  apparent  from  the  source  terra  in  Eq,  (II. 9),  which  depends  on  the 
definition  of  the  response  R.  This  highlights  the  complementary  nature  of 
forward  and  adjoint  methods.  The  forward  approach  gives  the  sensitivity  of 
all  responses  to  one  parameter  in  one  run,  whereas  the  adjoint  method  gives 
the  sensitivities  of  one  response  to  all  parameters  in  one  run.  From  a 
sensitivity  point  of  view,  the  adjoint  method  is  thus  clearly  more  applicable 
to  problems  with  large  numbers  of  input  parameters  and  few  responses- 
of-interest. 

5.4.3  Application 

A particularly  striking  example  of  the  power  of  adjoint  theory  in  the 
measurement  of  model  sensitivity  is  given  by  Oblow  [34]  . In  this  paper  Oblow 
applies  adjoint  theory  to  a linear  program  problem  and  to  other  nonlinear 
problems  with  and  without  given  constraints.  Oblow  concludes  that  "the 
developments  presented  make  it  clear  that  [adjoint]  sensitivity  theory  can  be 
extended  successfully  to  cover  a wide  class  of  algebraic  nonlinear  equations 
with  and  without  constraints."  He  continues,  "the  sensitivity  coefficients 
made  available  by  the  methods  developed  here  can  be  put  to  a number  of  uses. 
For  example,  Taylor  series  expansions  using  sensitivity  coefficients  can  be 
used  as  the  basis  for  a second-order  accurate  perturbation  theory  for  the 
nonlinear  systems  under  investigation.  In  addition,  a statistical  uncertainty 
analysis  of  system  responses  can  also  be  made  if  perturbation  results  are 
combined  with  assumptions  about  the  nature  of  the  uncertainties  in  the  system 
input  parameters". 
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5.5  Linear  Programming 


5.5.1  Description 

5. 5. 1.1  Postoptimality  Evaluation  The  method  of  linear  programming  is  so 
well  known  that  this  section  will  discuss  its  contribution  to  sensitivity 
analysis  methodology  with  only  a reminder  that  linear  programming  entails  the 
optimal  solution  of  an  algebraic  formulation  of  a problem  expressed  as  a 
linear  objective  function  subject  to  linear  constraints.  Its  connection  to 
sensitivity  analysis  is  the  postoptimality  evaluation  of  a solved  linear 
program. 

Following  Taha  [39],  the  changes  in  the  problem  include  variations  of: 

(1)  the  right-hand  side  of  the  constraints  (resource  vector),  (2)  the 
coefficients  of  the  objective  function  (cost  or  profit  values),  (3)  the 
coefficients  of  decision  variables  (matrix  of  technological  coefficients);  or 
the  addition  of:  (5)  new  variables,  and  (6)  new  or  secondary  constraints. 

The  sensitivity  associated  with  these  changes  is  determined  by  the  affect  they 
have  on  the  basic  variables  of  the  optimal  solution:  the  basic  variables  may 

remain  essentially  unchanged;  the  basic  variables  remain  but  their  values 
change;  or  the  basic  solution  changes  completely.  The  extent  to  which  these 
variations  can  be  made  without  changing  the  basic  variables  of  the  solution, 
or  essentially  changing  the  value  of  the  objective  function,  constitutes  a 
sensitivity  analysis  of  the  coefficients  and  constraints  of  a linear  program 
model.  A statement  of  the  primal  and  dual  formulations  of  a linear  program 
are  given  in  Appendix  B. 

5. 5. 1.2  Parametric  Sensitivity  A related  form  of  sensitivity  analysis  is 
called  parametric  linear  programming.  It  investigates  the  response  of  the 
solution  to  predetermined  linear  variations  in  the  coefficients  and  resource 
vector  based  on  the  primal-dual  relationships  used  in  the  postoptimality 
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formulation  [see  [39],  chapter  9].  Further  details  on  parametric  programming 
are  given  in  Appendix  B.  There  the  various  forms  an  analysis  may  take  are 
presented  and  a general  review  of  the  analytic  procedure  is  provided. 
Parametric  linear  programming  is  a well  known  technique  of  sensitivity 
analysis  and  has  been  incorporated  into  linear  programming  software  packages 
as  a standard  operating  procedure  for  linear  program  problems. 

5. 5. 1.3  Tolerance  Sensitivity  Of  considerable  interest  and  as  a innovative 

contribution  to  linear  programming  sensitivity  analysis  is  the  recent  work  of 
Wendell  [40]  on  the  subject  of  tolerance  evaluation  in  a linear  program 
problem.  His  abstract  summarizes  his  contribution:  "In  contrast  to 

'ordinary'  sensitivity  analysis  in  linear  programming,  the  tolerance  approach 
considers  simultaneous  and  independent  changes  in  the  objective  function 
coefficients  and  in  the  right-hand-side  terms.  This  approach  calculates  a 
maximum  tolerance  percentage  such  that  as  long  as  selected  coefficients  or 
terms  are  accurate  to  within  that  percentage  of  their  estimated  values,  the 
same  basis  is  optimal.  In  particular,  if  the  objective  function  coefficients 
are  accurate  to  within  the  maximum  tolerance  percentage  of  their  specified 
values,  then  the  same  solution  is  optimal."  Wendell  states  in  his 
introduction  that  his  tolerance  approach  is  a new  perspective  that  permits  a 
sensitivity  analysis  which  deals  with  more  than  one  coefficient  term  at  a time 
and  yields  the  largest  percentage  in  which  terms  may  vary  simultaneously  and 
independently  from  their  numerical  values  while  retaining  the  same  optimal 
basis.  His  method  also  yields  an  analogous  result  for  the  coefficients  of 
the  objective  function. 

5. 5. 1.4  Lagrangian  Parametric  Sensitivity  A method  of  considerable 
importance  for  the  study  of  parametric  sensitivity  has  been  developed  from  the 
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classical  Lagrangian  optimization  technique.  In  the  economic  literature  the 


method  is  referred  to  as  the  envelope  theorem,  but  its  derivation  and  use  in 
general  analytic  form  makes  it  universally  applicable  to  any  constrained  or 
unconstrained  objective  function.  The  theorem  provides  an  equation  for  the 
computation  of  the  rate  of  variation  of  an  objective  function  with  respect  to 
variation  in  a specified  parameter  of  interest.  The  following  derivation 
follows  Silberberg  [54]. 

Maximize  an  unconstrained  function  U,  defined  as: 

U - f ( xj  , • • . ( Xjjj  j a ) . 

From  the  Lagrangian  technique  for  the  optimum  value  of  U,  U* , we  write 

U*  = <J> (cx ) = f(xj*(a),  ...,  xn*(a)),  where  xj*(a)  are  the  optimal  values 

of  X£  as  functions  of  the  parameter  a.  The  sensitivity  coefficient  is, 

therefore: 

8U*  = = 8_f  • 

8a  8a  8a 

If  now  the  same  problem  is  solved  using  the  Lagrangian  technique  under  the 
constraint  functions  gj , i.e.: 

gj ( xj , ...,  xn;a)  = 0,  for  j=l , ...,  m,  then 
The  sensitivity  coefficient  becomes: 


Differentiating  each  constraint  function  with  respect  to  the  parameter  a,  we 
obtain: 


n 


8U*  = 8^  = J _8_f 


f 9xi*  + 8 f 


8 a 8a 

i=l 


n 


8x^*  8a  8a 

i=l 


£ 8gj  ^xi*  + ^ Sj  = 0,  for  each  j=l 


m; 
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3 <j> 

then  multiplying  by  Xj  and  adding  the  equations  to  3a  . 

The  equation  of  the  envelope  theorem  follows,  and  is  given  by: 

m 

3 u*  = 3L  = 3_f  + V“  Xj  _^gj_  , (5.5) 

3a  3a  3a  £_ 3a 

j = l 

m 

where  L is  the  Lagrangian  function,  f + \ Xjgj . 

j 

This  method  for  computing  the  sensitivity  coefficients  for  objective 
functions  suffers  from  the  same  computational  difficulties  associated  with  the 
Lagrangian  optimization  technique,  primarily  because  the  optima]  values  and 
the  multipliers  may  be  tedious  to  obtain  for  complex  or  large  problem 
formulations . 

However,  the  method  is  attractive  because  it  computes  the  sensitivity  of 
objective  functions  with  respect  to  any  parameter  in  the  formulation;  that  is 
the  parameter  may  be  a resource  level,  a cost  or  profit  coefficient,  or  a 
technical  coefficient.  Since  the  objective  function  and  constraint  functions 
may  be  nonlinear  the  envelope  theorem  permits  the  concept  of  linear  parametric 
sensitivity  analysis  to  be  carried  over  into  the  the  domain  of  the  parametric 
sensitivity  of  nonlinear,  optimization  problems. 

5.6  Sensitivity  Functions 
5.6.1  Description 

The  sensitivity  function  is  derived  from  a differential  equation 

F(x,x,x,t ,q)=0  (5.6) 

in  which  q is  a parameter  to  be  perturbed.  If  we  assume  a solution  to 
equation  (5.6),  x=x(t,q),  we  next  determine  the  effect  on  x of  a 
perturbation,  Aq,  in  q.  This  gives 

F ( x , x , x , t , q+A  q ) =0  (5.7) 
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as  the  new  system  differential  equation.  We  now  define  a variable,  termed  the 
sensitivity  coefficient,  as  u(t,q)  which  is  the  derivative  of  x with  respect 
to  q : 

u(t,q)  = lira  x(t ,q+Aq)— x(t ,q)  _ dx(t,q)  (5.8) 

Aq+0  Aq  dq 

The  sensitivity  coefficient  measures  the  change  in  x due  to  a 
pertrubation  in  q.  Differentiating  equation  (5.6)  with  respect  to  q gives: 


3 F 

9x  + 

3F  3x 

+ 3F  3x 

+ 

3 F = 

Tx 

3q 

31  3q 

3x  3 q 

3 q 

with  3x 

= 3 

9x  = u 

rsC* 

and  ii2S  = 

32 

3x  = 

3q 

3 1 

3q 

3q 

Tt2 

3 q 

then  equation  (5.8)  takes  the  form 


(5.9) 


l£ff  + iF5  + 3Fu  + 3F  = 0 • (5.10) 

3x  3x  3x  3q 

Equation  (5.11)  is  the  sensitivity  equation  which  is  to  be  solved  for  u(t,q). 
The  sensitivity  equation  is  a linear  equation,  and  provided  the  initial 
conditions  for  equations  (5.6)  and  (5.7)  are  the  same,  all  initial  conditions 
of  equation  (5.11)  are  equal  to  zero.  See  References  [41]  and  [42]  for  more 
detail  on  the  theory  and  scope  of  the  sensitivity  function.  Reference  [43] 
demonstrates  an  application  to  a world  dynamics  model  based  on  a system  of 
nonlinear  equations  of  the  first  order. 
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5.7  Eigenvalue  Functions 
5.7.1  Description 

In  a system  of  linear  first-order  differential  equations 
y = Ay  (5.12) 

where  y is  an  n-component  state  vector  and  A an  n by  n matrix  of  coefficients, 
the  solution  is  of  the  form 

y = cxe^t  (5.13) 

where  x is  an  n-component  vector  and  c is  an  arbitrary  scalar.  When  equation 
(5.13)  is  substituted  into  equation  (5.12)  the  problem  is  reduced  to  finding 
the  eigenvalues  of  the  system  of  equations. 

Ax  = A x (5.14) 

If  it  is  necessary  to  study  the  parametric  effects  of  variations  in  the 
eigenvalues  then  equation  (5.14)  is  partially  differentiated  with  respect  to  a 
parameter  q and  the  variational  equation  is  obtained 

3A  x±  +A^fi  = Ai^i+  9xi  Xj  (5.15) 

3 q 3 q 3 q 3 q 

for  each  i,  the  index  over  the  vector  x.  After  some  mathematical  manipulation 
the  sensitivity  coefficient  is  expressed  by 

(*±  xi>  vi^ 

^_i  / (5.16) 

3 q (xi}  Vi) 

where  the  notation  (*,*)  indicates  a scalar  product,  and  Vi  is  the  ith 
eigenvector  of  the  system  when  the  matrix  A is  transposed.  For  a given 
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eigenvalue,  X^,  the  eigenvectors  and  v-^  are  calculated,  then  equation 
(5.15)  is  used  to  compute  the  sensitivity  coefficient.  This  description  is 
based  on  material  contained  in  Reference  [42]. 

5.8  Large  Error  Sensitivity  Analysis 
5.8.1  Description 

By  large  error  we  mean  those  errors  of  sufficient  magnitude  which 
render  the  classical  linear  error  equations  invalid  as  measures  of 
sensitivity.  Classical  error  analysis  assumes  that  variations  behave 
with  necessary  properties,  (1)  that  they  are  small,  i.e.,  neighborhood 
variations,  and  (2)  that  on  average  there  are  as  many  negative  variations  as 
there  are  positive  variations.  This  second  assumption  is  critical  to  the 
classical  formulation  of  variance. 

The  evaluation  of  sensitivity  coefficients  for  large  errors  rests  on 
using  more  complicated  formulae  where  second  order  terms  are  not  dropped,  on 
using  the  variational  form  of  the  differential  method  described  in  Section 
3.3,  or  on  the  use  of  interval  analysis  in  which  extreme  values  of  a function 
serve  to  delineate  the  range  of  a function  under  perturbation. 

Rahman  [44]  gives  several  examples  of  simple  operations  such  as  x/y  where 
both  variables  are  subject  to  variation.  Letting  z = x/y,  the  largest  change 
in  x/y,  excluding  y=0,  will  be  when  x is  at  its  largest  value  and  y at  its 
smallest.  Similarly  the  smallest  value  of  Z will  occur  when  x is  minimized 
and  y takes  its  greatest  variational  value.  This  gives  Zu  and  Z l where 

Zu  = min  max  [Z=x/y] 
y x 

and  Z^  = max  ®in  [Z=x/y] 
y x 
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or  Zl  < Z < Zu.  If  anything  is  known  about  the  distribution  of  x and  y 
in  the  intervals  of  their  variation  then  the  theory  described  in  Section  4.3 
can  be  used  to  obtain  estimates  of  statistical  parameters  for  Z. 

This  concept  can  be  generalized  by  estimating  the  absolute  maximum  and 
minimum  values  of  a function  over  the  permitted  hypercube  of  variations  of  the 
independent  variables.  If  the  range  of  variation  is  large  but  the  relative 
error  (ie.  6x/x)  is  less  than  one  then  some  simplifying  substitutions  may  be 
possible  in  the  determination  of  the  maximum  and  minimum.  If  the  function 
does  not  possess  a workable  probability  distribution  over  the  hypercube  of 
permitted  variation  then  the  sensitivity  can  be  expressed  as  an  absolute  error 
(±  e),  or  as  a percent  error  of  this  expected  value  of  the  function  over  the 
range  of  the  dependent  variable. 
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6.0  SPECIAL  ANALYTICAL  METHODS  OF  SENSITIVITY  ANALYSIS 


6.1  Introduction 

In  this  section  three  very  specialized  methods  of  sensitivity  analysis 
are  described.  They  are  included  for  various  reasons  - the  Fourier 
coefficients  method  has  some  interesting  properties,  geometric  programming  is 
a useful  tool  for  studying  engineering  costs  prior  to  project  commitment,  and 
Catastrophe  Theory,  much  like  Game  Theory,  provides  a conceptual  structure  in 
which  many  problems  of  a diverse  nature  can  be  analyzed  with  greater  insight. 

6.2  Fourier  Coefficients 
6.2.1  Description 

The  principal  source  is  Cukier,  Levine,  and  Shuler  [45]  from  which  we 
abstract,  with  minor  editing,  a portion  for  the  description  of  the  method. 

"Separate  in  definition  from  the  independent  variables  and  parameters  are 
the  fixed  constants  of  a model,  which  do  not  vary  within  the  context  of  the 
class  of  problems  of  interest  to  the  model  user,  and  whose  values  can  be 
precisely  specified.  It  should  be  noted,  of  course,  that  what  is  a fixed 
constant  in  the  context  of  one  situation  might  be  a parameter  in  the  context 
of  another  situation;  the  distinction  depends  on  the  particular  case  on  hand. 
The  fact  that  the  parameters  can  take  on  a range  of  values  suggests  that  a 
statistical  approach  to  sensitivity  analysis  is  appropriate.  Instead  of 
considering  the  effect  on  the  output  functions  of  one-at-a-time  variations  in 
each  of  the  parameters,  as  in  a "brute  force"  method,  we  will  construct 
outputs  averaged  in  one  operation  over  probability  distributions  of  all  the 
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parameters.  The  distribution  of  the  parameters  can  arise  because  of 
experimental  uncertainties  or  theoretical  approximations,  because  of 
"ignorance"  of  the  value  within  certain  reasonable  bounds,  or  might  represent 
upper  and  lower  limits  due  to  "stops"  on  the  physical  controls  of  the  systems 
being  modeled. 

"Our  method  of  sensitivity  analysis  proceeds  by  relating  the  probability 
distribution  of  each  parameter  to  a frequency  and  one  new  parameter  s which, 
as  s varies,  carries  all  the  parameters  through  their  range  of  variation.  The 
parameter  s is  varied,  and  the  frequencies  are  chosen  in  such  a way  that  the 
output  variables  at  any  fixed  time  become  periodic  in  s. 

"The  output  variables  can  then  be  evaluated  using  Fourier  Analysis.  As 
we  shall  show  below,  the  Fourier  coefficients  represent  an  average  of  the 
output  variables  over  the  uncertainties  of  all  the  parameters.  A unique 
correspondence  between  the  Fourier  coefficients  for  the  frequency  w^  and  all 
its  harmonics  and  the  sensitivity  of  the  output  variables  to  the  kth  parameter 
is  established.  We  compress  all  this  information  into  partial  variances  Sw^ 
which  are  the  normalized  sums  of  the  squares  of  the  Fourier  coefficients  of 
the  fundamental  frequency  w^  and  all  its  harmonics.  If  Sw^  < Sw£  for  a given 
output  variable,  then  this  output  variable  is  less  sensitive  to  the  kth 
parameter  than  to  the  jth  parameter.  Thus,  the  partial  variances  measure  the 
average  sensitivity  of  an  output  function  to  the  variation  (or  uncertainty)  of 
a particular  parameter.  This  average  is  over  the  range  of  uncertainties  of 
all  the  parameters,  with  their  appropriate  probability  distributions,  with  the 
exception  of  the  parameter  being  considered.  For  this  parameter,  the 
statistical  property  calculated  is  the  variance. 
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"The  sensitivity  analysis  is  nonlinear  so  that  it  permits  us  to  examine 
large  deviations  from  the  nominal  parameter  values.  In  addition,  since  all 
parameters  are  varied  simultaneously,  one  explores  regions  of  parameter  space 
where  more  than  one  parameter  is  far  from  its  nominal  value.  Because  of  this 
thorough  and  systematic  exploration  of  the  parameter  space,  it  often  turns  out 
that  sensitivities  of  an  unexpected  nature  are  revealed.  A careful  study  of 
the  model  will  then  reveal  some  complex  coupling  between  variables,  unexpected 
prior  to  the  analysis,  which  leads  to  the  observed  sensitivity.  In  this 
fashion,  one  obtains  deeper  insights  into  the  structure  of  the  complex  system 
being  studied.  Another  frequent  and  important  finding  is  that  a number  of 
senstivity  coefficients  corresponding  to  a large  set  of  parameters  turn  out  to 
be  negligible.  This  permits  one  to  reduce  the  complexity  of  the  set  of  model 
equations  and  focus  one's  attention  on  a greatly  reduced  set  of  equations." 

6.3  Geometric  Programming 
6.3.1  Description 

Geometric  programming,  a well  known  method  of  optimization  of  nonlinear 
problems,  gives  the  optimal  cost  or  profit  before  a corresponding  design  or 
plan  is  implemented.  The  method  is  described  in  detail  in  Duff in,  Peterson 
and  Zener  [46]  , and  is  given  as  excellent  exposition  in  Wilde  and  Beightler 
[47]. 

Instead  of  obtaining  an  optimal  solution  of  the  decision  variables 
initially,  geometric  programming  finds  the  optimal  distribution  of  the  total 
cost  among  the  terms  of  the  objective  function,  and  when  these  optimal 
assignments  are  obtained  the  optimal  cost  follows  directly  from  a 
computation.  Once  the  optimized  cost  and  weights  are  determined,  they  can  be 
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used  to  determine  the  optimal  policy.  Geometric  programming  can  be  used  to 
study  the  policy  consequences  of  variations  in  cost  and  in  engineering  design 
variables. 

6.4  Catastrophe  Theory 
6.4.1  Description 

Catastrophe  theory  describes  the  phenomena  associated  with  sudden  changes 
of  state  in  a system  caused  by  smooth  alterations  in  the  control  parameters  of 
the  system.  The  theory  is  not  easy  to  master.  It  is  derived  from  topology 
which  is  needed  to  describe  forces  as  smooth  surfaces  of  equilibrium,  but  it 
also  is  not,  according  to  Poston  and  Steward  [48],  a single  thread  of  ideas, 
but  embraces  geometry,  algebra,  singularity  theory,  physical  intuition  and 
experiment.  The  sources  on  this  subject  are  References  [48-52], 

The  underlying  concept  of  catastrophe  theory  is  the  recognition  of  the 
various  kinds  (in  a mathematical  sense)  of  catastrophes  and  the  relationship 
between  the  catastrophe  manifold  and  the  control  parameter  space. 

The  catastrophe  manifold  or  equilibrium  surface  is  a set  of  points 
(x,{ai)^),  where  x is  a state  variable  and  {a^j^  is  a set  of  N control 
variables,  a^.  The  map  of  (X,{a-£}ig)  will  have  one  or  more  folds  in  its 
surface  which  are  a function  of  the  potential  of  the  energy  of  the  system. 

The  projection  of  the  equilibrium  surface,  {a^}j^  or 

(x,{ai}N)  > {a-^  (5.17) 

defines  the  control  parameter  space.  The  folds  of  the  equilibrium  surface  map 
into  cusps  in  the  control  parameter  space.  As  long  as  the  values  of  { a-^ > are 
outside  the  enclosed  region  of  the  cusps  the  state  variable  will  experience 
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smooth,  stable  transitions.  When,  however,  the  values  of  {a^}  cross  into  the 
cusp  region  the  state  variable  is  immediately  unstable  and  will  jump  to  a 
higher  or  lower  value  when  the  {a^}  path  or  trajectory  exits  from  the  cusp 
region.  Throughout  this  process  the  path  of  the  control  parameters  values  is 
smooth  and  continuous. 

This  subject  bears  on  both  the  qualitative  and  quantitative  nature  of 
sensitivity  analysis.  The  qualitative  aspect  lies  in  the  insistence  the 
theory  makes  on  the  model  builder  to  study  his  subject  closely  for  unstable 
or  jump  conditions  in  the  state  variables  of  his  system,  and  to  include  these 
possibilities  in  his  model.  Far  too  often  models  fail  because  the  underlying 
process  was  assumed  to  be  a steady-state  system  when  in  fact  it  was  subject  to 
"catastrophic”  changes  of  state. 

The  quantitative  aspect  of  catastrophe  theory  lies  in  its  broad 
application  to  many  diverse  types  of  problems.  It  has  been  used  in  physics, 
biology,  social  behavior,  species  population  analysis,  medical  problems  and 
models  of  ecological  phenomena.  It  is  not  necessary  in  all  cases  to  use 
catastrophe  theory  to  assess  systemic  instability.  However,  even  in  these 
cases  the  knowledge  of  its  elements  should  be  very  helpful  in  understanding 
the  phenomena  of  state  transitions.  There  are  several  catastrophe  "flags" 
which  signal  certain  conditions.  Some  aspects  of  system  dynamics  and 
associated  flags  are  discussed  in  Appendix  C.  An  understanding  of  systemic 
flags  enhances  the  quality  of  a model  if  the  model  builder  anticipates  the 
impact  of  these  situations  and  incorporates  their  consequences  in  the  model 
structure  and  dynamics. 

Catastrophe  theory  can  contribute  to  a quantitative  sensitivity  analysis 
by  using  its  basic  ideas  to  compare  model  results  with  historical  results. 
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This  application  is  particularly  relevant  to  forecast  models,  economic  models 
and  time-series  analysis.  A plot  of  the  values  of  the  control  variables  over 
time  versus  system  state  data  will  show  the  effects  of  the  model  performance 
in  comparison  to  the  system  performance  over  the  period  of  interest.  This 
comparison  should  provide  a revelation  of  detail  that  will  indicate  how  and  in 
what  manner  control  parameters  relate  to  the  system  state  variables  and 
to  the  simulation  of  the  system. 
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7.0  MODELING  AND  SENSITIVITY  ANALYSIS 


7.1  Introduction 

In  this  final  section  various  ideas  will  be  detailed  which  bear  on  the 
need  to  bring  sensitivity  analysis  into  a closer  and  more  substantive 
relationship  with  model  development.  Since  most  of  the  methods  covered  in 
this  survey  are  external  to  the  model  they  result  in  a large  investment  in 
time,  in  data  preparation,  in  input  selection,  and  in  runs  of  the  model.  Some 
of  these  separately,  and  undoubtedly  all  of  them  together,  either  defeat  or 
discourage  sensitivity  studies  because  of  the  burden  of  having  to  impose  the 
sensitivity  methodology  on  the  model  in  an  external  and  inefficient 
procedure. 

Attempts  will  be  made  to  suggest  ways  to  alleviate  some  of  these  burdens 
of  testing  a model's  sensitivity  to  input  data  uncertainty  and  control 
parameter  variation. 

7.2  Modeling  Control  Parameters 

A parameter  is  a value  of  a constant  that  determines  the  character  of 

behavior  of  something  which  is  a function  of  that  parameter.  The  value  of  the 

control  parameter  is  determined  at  the  discretion  of  the  user  of  the  model. 

Most  models  have  many  control  parameters  embedded  in  their  logic  which 
are  usually  constants  whose  values  are  set  initially  as  input  data  to  the 

model.  The  parameter  is  then  fixed  at  its  given  initial  value  for  the 

entirety  of  the  model  run.  Static  models,  for  the  most  part,  are  adequate 
which  use  this  procedure,  but  models  that  simulate  systems  that  are  dynamic  or 
predictive  are  severely  restricted  by  this  approach. 
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To  address  this  problem,  and  as  appropriate  for  the  model  in  question, 
each  important  control  parameter  should  be  constructed  as  a time  or  index- 
dependent  variable,  with  each  value  per  time-period  or  index  number  given  as 
input.  This  allows  the  decision  maker  to  vary  the  control  parameter  according 
to  the  conditions  to  be  examined  and  permits  a more  realistic  simulation  of 
exogenous  events  which  impact  on  or  influence  system  behavior. 

Another  version  of  this  idea  is  to  assign  values  to  a control  parameter 
vector  and  have  the  model  access  the  appropriate  value  of  the  control 
parameter  successively,  according  to  an  index  or  a prescribed,  sequential 
assignment. 

Making  the  control  parameter  a variable  by  one  of  the  above  techniques 
provides  the  user  with  a program  procedure  for  studying  dynamic  changes  in  the 
model  scenario.  It  also  provides  a method  for  perturbing  the  values  of 
control  variables,  in  a variational  sense,  that  will  generate  data  which  can 
be  used  in  a sensitivity  analysis  of  the  model's  performance. 

7.3  Modeling  Sensitivity  Methodology 

In  almost  every  scheme  of  current  methodology  for  analyzing  the 
sensitivity  of  a model  the  emphasis  is  either  on  appropriate  generation  of 
input  data,  defining  response  functions,  or  performing  a postpartum 
statistical  analysis  in  which  the  model  serves  as  the  experiment.  Other 
methods  analogize  the  model  as  an  operator  performing  on  input  data  to  produce 
output  which  is  then  tested  on  the  basis  of  linear  sensitivity  formulations. 
The  model  is,  in  these  experiments,  perceived  and  treated  as  a black  box.  The 
sensitivity  analyses  of  this  approach  provides  an  overall  measure  but  it  does 
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not  describe  the  relational  dependencies  and  sensitivities  experienced  within 
the  components  of  the  model,  of  which  there  may  be  many  and  of  a complex 
nature. 

In  order  to  reduce  some  of  these  objections,  and  at  only  a reasonable 
increase  in  the  initial  investment,  it  is  suggested  that  sensitivity  analysis 
should  be  an  integral  part  of  a model,  incorporated  into  its  logic,  and 
structured  to  provide  information  on  validation.  The  general  concept  of  this 
integration  would  be  as  follows. 

A model,  in  its  initial  stages  of  development,  would  normally  be  divided 
into  its  major  components,  as  most  models  are;  but  this  partitioning  would  be 
extended,  within  the  components,  to  groupings  of  self-contained  computations, 
which  we  will  call  mathematical  modules.  Assuming,  for  simplicity,  that  the 
flow  chart  of  a portion  of  the  model  is  shown  as  in  Figure  1. 


Figure  1 

It  is  not  important  from  this  point  on  to  retain  the  identity  of  the 
components  since  the  design  is  based  only  on  the  modules.  The  scheme  is  now 
enhanced  by  providing  each  module  with  a program  which  provides  a sensitivity 
subroutine  appropriate  for  the  computations  performed  by  the  module.  This 
enhances  the  structure,  as  shown  in  Figure  2. 
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Figure  2 

Next  it  is  necessary  to  provide  two  operational  paths  through  the  flow  chart 
of  Figure  2.  The  path  shown  in  Figure  3 is  the  normal  operational  use 
of  the  model  without  invoking  the  sensitivity  subroutines  (SS). 


To  obtain  the  second  path  it  is  necessary  to  construct  a control  vector,  "x" , 
composed  of  a sequence  of  logical  control  switches,  one  for  each  SS,  which 
which  will  be  "opened"  or  "closed"  according  to  whether  all  or  some  of  the 
sensitivity  subroutines  will  be  activated.  If  they  are  all  used  then  the 
operational  flow  will  be  as  in  Figure  4. 


Figure  4 presents  a basis  for  a number  of  variations  on  the  essential 
concept.  For  example,  if  each  senitivity  subroutine  can  be  optioned  to  output 
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its  calculations,  then  a running  profile  of  the  modular  sensitivities  can  be 
obtained.  The  sensitivity  subroutines  could  be  written  to  pass  their 
computations  forward  to  the  next  SS  operation,  and  so  accumulate  their  results 
for  some  appropriate  output  option  point.  The  information  needed  by  each  SS, 
primarily  the  variations,  could  be  stored  directly  in  an  SS  program,  or  some 
scheme  could  be  programmed  to  have  each  activated  SS  compute  its  own 
variational  values  based  on  statistical  parameters,  relative  errors,  or 
functions . 

The  control  variable  X would  be  a vector  of  "on"  or  "off"  signals  which 
would  be  processed  as  the  initial  operation  of  the  model.  This  vector  would 
determine  the  process  path  of  the  model  for  that  run.  Another  control  vector 
could  be  used  to  signal  to  each  SS  how  it  is  to  compute  the  values  of 
the  variations  it  is  to  use  in  computing  the  sensitivity  response  of  its 
associated  mathematical  module. 

By  whatever  means  it  is  accomplished,  by  the  above  concept  or  by 
inventions  of  new  techniques,  it  appears  that  with  the  growing  interest  in 
model  validation  and  sensitivity,  it  is  essential  to  incorporate  sensitivity 
methodology  directly  into  the  model  structure. 

The  current  practices  of  performing  sensitivity  analyses  are  not  wholly 
satisfactory  because  of  the  assumptions  which  must  be  made  and  the  often  high 
cost  of  using  a model  as  a black  box  in  an  experimental  design  from  which  the 
sensitivity  analysis  is  carried  out  as  a postpartum  exercise.  Incorporating 
the  sensitivity  analysis  into  the  model  logic  provides  more  detail  to 
validation  and  sensitivities  studies,  and  offers  the  user  a procedure  for 
assessing  the  structural  details  and  dynamics  of  the  system  and  its 
simulation. 
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Appendix  A 


Part  I:  Summary  of  Principal  Error  Formulas. 

The  principal  formulas  associated  with  the  propagation  of  errors  are 
presented  for  convenience  and  for  the  purposes  of  handy  reference.  These 
formulas  are  given  in  functional  form  in  order  to  show  their  character  and 
structure,  and  to  generalize  the  formulas  to  multi-dimensional 
representations . 

In  the  formulas  the  following  definitions,  symbols,  and 
notations  will  be  used: 

6x  is  used  to  indicate  a variation  in  a prescribed  variable,  here 
designated  as  x; 

F is  used  to  indicate  a function;  F(x^)  is  used  to  indicate  a function  of 
a set  of  variables  {x^},  with  the  index  i running  from  1 to  k; 

rx^xs  is  the  correlation  coefficient  of  variables  x^  and  xs; 

n is  used  to  designate  the  size  of  the  sample; 

c?x2  is  the  variance  of  variable  x^; 


i 


xA  is  the  average  of  the  variable  x^; 

°x^xs  is  the  covariance  between  x^  and  xs; 


i,  j,  s are  dummy  indices  used  in  summation  terms 


Absolute  Error 


The  absolute  error  of  a function  W=F(x^)  is  defined  as 


( A — 1 ) 


where  8F  are  evaluated  from  some  given  point 
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Relative  Error 


The  relative  error  of  a function  W=F(xi)  is  defined  as 


6w  = I [ 3FVxi  , 

w i=i  \ ax±;  w 
Probable  Error 

The  probable  error  of  a function,  given  an  error  e„  in 

i 

each  x^,  is  defined  as: 


(A-2) 


the  measurement  of 


Mean  of  a Function 

The  mean  of  a function  W=F(xi)  is  given  as: 
W = F(xt)  + R 


(A-3) 


(A-4) 


Where  R is  equal  to: 


R 


k-1 

+ l 


j-i 


a 

xixs  . 
n 


(A-5) 


When  n becomes  large  R tends  to  zero,  leaving 

W = F(Xi)  (A-6) 

as  the  general,  approximating  expression  for  the  function  mean  where  the 
variations  are  assumed  small  and  normally  distributed.  No  assumption  of 
independence  is  made  in  (A-5);  if  independence  pertains  then  (A-5)  simplifies 
to 
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(A-7  ) 


R 


= 1 I 


a2F 


2 i=l  l 3x| 


JiL 

n 


Variance  of  a Function 

The  variance  of  a function  w=F(x^)  is  defined  by  the  equation 


k 

l 

i=1  \ ax. 


°F2  = l °xi2  + 2 l l 1 rXj  xs  axj°! 


j = l s=j+l  \ 3X . 3X 
J s 


(A-8) 


where  ax^  and  °xs  are  standard  deviations. 

Covariance  of  Functions 

Given  that  U=F(x£)  and  V=G(x-j[)  are  two  functions  their  covariance  is 
defined  as: 


°uv  = l 
i=l 


au  3v 


3x.  3x. 

i i 


k-1 


,xi  + I l I1.  + 1L  —Yx jxs0xj0xs  ( a-9) 


j = l s=j+l  \ 9x>  9X 
J s 


3x  3x. 
s J 


Systematic  Error 

A systematic  error  is  defined  as  a fixed  deviation  which  is  in  each 
measurement  of  a variable  in  a particular  sequence  of  measurements.  Given 
W=F(x£),  the  systematic  error  is: 


6W|  = I | 9F  SxJ 


i=l  9xi 


(A-10) 


The  user  of  these  formulas  should  keep  in  mind  that  they  are 
approximations,  and  in  strict  mathematical  terras  the  equality  sign  is  not 
literally  true.  Further,  the  underlying  assumptions  on  which  these  formulas 
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are  based  should  be  reviewed  for  compliance  to  statistical  requirements  in  any 
problem  of  sensitivity  analysis  in  which  they  are  used. 

The  in  the  above  equations  are  not  restricted  to  variables;  they  can 
be  thought  of  as  parameters,  coefficients,  or  input  data  for  which  a 
sensitivity  response  is  to  be  measured. 
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Part  II:  Statistics  of  the  Index  w=x/y. 

The  theory  underlying  the  statistics  of  w=x/y  is  presented  in  a thorough 
and  exhaustive  analysis  by  Fieller  [6]  which  is  far  too  difficult  and  lengthy 
to  reproduce  here.  The  general  idea,  based  on  [6] , will  be  reviewed  for  those 
interested  in  the  background  of  these  approximations. 

Following  Fieller,  the  ratio  x/y  is  expanded  by  their  variations,  fix  and 

fiy : 

Let  x = x + fix  and  y = y + fiy,  x is  the  expected  value  of  x;  y is  the 
expected  value  of  y;  expand  as  follows: 


U - x+6x  _ x 

(i+sA 

i-  Ax 

+ Azi 

- Axi  + ... 

y4-5y  “ 

-) 

y 

y 

■ J 

Then  compute  the  n-th  moment  about  zero  of  the  distribution  of  w,  from  the 

! 

equation  for  W1.  To  obtain  this  value  (based  on  a derivation  by  Merrill,  see 
reference  [6] , page  436)  retain  the  products  of  (fix)r  (6y)s  as  far  as  the 
eighth  order,  and  take  for  their  mean  values  the  product-moments  of  the  normal 
surface.  This  holds  only  for  |6y|  < y,  but  it  is  valid  when  applied  to  the 
interior  of  the  probability  contour.  Merrill's  values  of  the  moments  may  be 
taken  as  the  moments  of  the  distribution  of  w in  a curtailed  normal 
population.  This  provides,  then,  the  basis  for  the  formulas  given  in  Section 
3. 3. 1.1,  eauaitons  (3.2)  and  (3.3). 

;i 
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Appendix  B 


Parametric  Programming:  Sensitivity  Analysis 

Parametric  linear  programming  can  involve  the  evaluation  of  six  types  of 
modifications.  Given  a typical  linear  program  as: 
x > 0 y > 0 


Ax  < d A’y  > c 

max  c'x  = z min  d'y  = z 

(Primal)  (Dual) 

In  these  formulations  A is  the  matrix  of  technical  coefficients,  d is  the 
vector  of  resources,  c is  the  vector  of  costs  or  profits,  x is  the  vector  of 
allocations,  and  z is  the  objective  function,  which  is  a scalar  product.  The 
symbol  indicates  the  transpose  of  the  matrix  A or  the  vectors  c and  d. 

Two  properties  which  link  the  primal-dual  solutions  together  are  important  to 
parametric  sensitivity  analysis.  The  first  property  states  that  the  solution 
of  the  dual  problem  acts  as  an  upper  bound  on  the  solution  of  the  primal 
problem.  The  second  property  states  that  if  B is  the  basic  matrix  of  the 
solution  to  the  primal  problem  then  the  optimum  solution  of  the  dual  is 
Y = cB“l , where  c is  the  cost  vector  associated  with  the  primal  solution. 

The  impact  of  a parametric  sensitivity  analysis  on  the  integrity  of  these 
two  properties  should  be  assessed  with  regard  to  the  consistency  the 
modifications  have  with  the  primal-dual  relationship.  A parametric 
sensitivity  analysis  of  a linear  program  may  take  any  of  the  following 
modifications . 
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The  six  possible  modifications  are: 

(1)  c changes  by  a discrete  amount, 

(2)  d changes  by  a discrete  amount, 

(3)  a single  column  of  matrix  A changes  by  some  amount, 

(4)  a single  row  of  A changes  by  some  amount, 

(5)  a constant  is  added  to  the  system,  and 

(6)  a variable  is  added  to  the  system. 

These  six  cases  are  treated  thoroughly  by  Simonnard^  and  the  reader  is 
referred  to  this  text  for  an  excellent  exposition  and  analysis  of  these 
subjects . 

Only  the  first  modification,  changes  in  the  cost  vector  c,  will  be 
studied  in  this  Appendix.  The  analysis  entails  calculating  the  bounds  in 
which  c may  vary  without  altering  the  optimal  solution.  If  the  values  of  c, 
say  c^,  are  varied  by  a constant  value  t then  the  objective  function  in  the 
primal  becomes  (c'+t)x.  The  parameter  t will  appear  in  the  tableau  in  the 

^Simonnard,  M. , Linear  Programming,  Prentice-Hall,  Inc.,  1966,  translated  by 
William  S.  Jewell,  Chapter  7. 
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last  row  and  the  simplex  method  will  yield  several  inequalities  on  t. 

These  inequalities  define  the  interval  on  t for  which  the  optimal  solution 
will  not  be  altered.  The  objective  function  is  then  defined  as 
(c'+  tj^)  x*  < z < (c'+tu)x*,  where  x*  is  the  optimal  solution.  If  t is 
outside  the  interval  [t^,  tu]  then  the  optimal  solution  is  changed. 

If  the  problem  is  to  study  variations  in  the  requirements  vector  d,  the 
linear  program  may  be  recast  in  the  dual  form  and  variations  in  d may  be 
examined  in  the  same  manner  as  was  done  for  c in  the  explanation  above. 

Kreko*  develops  a general  theory  for  variations  in  c or  d.  He  defines  t^ 
and  tu  as  characteristic  points  of  the  parameter  t and  elaborates  the 
conditions  for  which  optimal  solutions  change  as  a function  t.  His  analysis 
describes  the  solutions  associated  with  the  characteristic  points  of  the  t 
domain,  and  develops  a procedure  for  establishing  upper  and  lower  bounds  on  t 
and  the  associated  optimal  solutions  associated  with  a monotonic  sequence  of 
characteristic  points. 


^Kreko,  B. , Linear  Programming,  American  Elsevier  Co.,  Inc.,  1968,  Sections 
4 . 6 and  12.3. 
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Appendix  C 


Catastrophe  Flags  and  System  Dynamics 

The  text  of  Section  6.4,  Catastrophe  Theory,  is  augmented  here  with 
greater  detail  on  the  dynamics  of  system  stability  as  a function  of  changes  in 
control  parameters,  and  with  a discussion  of  the  catastrophe  flags,  which  have 
a bearing  on  model  dynamics,  design  and  validation. 

Figure  C-l  depicts  the  geometry  of  the  cusp  catastrophe  associated  with  a 
parabolic  force  potential  manifold.  The  response  surface,  M,  is  the  folded 
sheet  in  three  dimensions.  Its  projection,  C,  is  the  control  parameter  space, 
which  is  shown  in  Figure  C-l  with  two  parameters  (a,b)  and  a cusp  emanating 
from  the  point  P.  Four  regions  in  the  control  space  are  marked  as  , R2 , R3 , 
and  R^.  Similarly  in  the  manifold  space  are  regions  designated  as  Q0,  Oj  , and 
Q 2.  The  points  u and  v in  M are  "jump"  points  at  which  the  state  variable 
jumps  to  only  stable  minimum  values.  Regions  and  O2  have  minimum  potential 
values.  The  region  of  the  fold  has  maximum  potential  value  and  is 
unstable . 

The  discussion  which  follows  on  the  catastrophe  flags  includes  the 
general  phenomena  of  system  dynamics  presented  in  Section  6.4. 

1 . Modality  Flag 

Modality  refers  to  the  condition  that  the  system  has,  in  the  case 
depicted  in  Figure  C-l,  two  or  more  distinct  physical  states  when  the  control 
parameters  have  values  lying  within  the  cusp-shaped  regions.  The  projection 
of  any  point  lying  within  the  cusp-line  boundary  will  pierce  the  manifold  in 
three  places.  All  such  points  are,  therefore,  states  of  instability  for  the 
system. 
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2.  Sudden  Jump  Flag 


A sudden  change  in  the  state  of  the  system  occurs  when  the  state  variable 
jumps  from  one  sheet  of  the  manifold  to  the  other.  In  Figure  0-1 , this  occurs 
when  the  control  parameters  describe  a path  from  regions  Rj  or  R2  to  R3  or  R4 
where  the  path  enters  the  cusp  region  and  exits  by  crossing  a cusp  line.  A 
sudden  jump  can  occur  when  the  control  parameters  move  the  state  variable  from 
Q2  and  Qi  by  traversing  through  the  point  v or  the  cusp-line  (0o,v).  Moving 
the  state  variable  from  region  to  O2  by  crossing  the  cusp-line  (Q0,  u) 
causes  a similar  sudden  jump  to  equilibrium  in  the  Q2  region.  A catastrophic 
jump  occurs  when  a smooth  and  continuous  variation  of  control  parameters 
causes  a discontinuous  change  of  state. 

3.  Divergence  Flag 

A divergence  is  said  to  occur  when  a smooth  and  continuous  change  in 
control  parameters  leads  to  a smooth  and  continuous  change  in  the  system  state 
variable.  This  can  be  seen  in  Figure  C-l.  If  the  system  state  variable  is  in 
the  region  near  Qq  it  can  be  moved  into  region  O2  when  the  parameters  follow  a 
smooth  path  from  Rj  to  R2.  A similar  movement  can  lead  the  state  variable 
from  Q0  to  Qj  as  the  control  parameters  move  along  a path  from  Rj  to  R3  (or 
R4)  which  does  not  traverse  the  cusp  region.  So  by  smooth  transitions  the 
state  variable  can  be  brought  to  divergent  values,  one  in  the  region  of  Q2 , 
the  other  in  the  region  Qj . 

4 . Hysteresis  Flag 

This  condition  is  included  because  it  may  have  some  bearing  on  the 
evaluation  and  behavior  of  modeling  design.  Hysteresis  describes  the 
situation  in  which  if  the  path  of  the  control  parameters  is  reversed  the  path 
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of  the  state  variable  is  not  necessarily  reversed.  Reversibility  is  possible 


if  the  parameter  path  does  not  exit  from  the  cusp  region.  The  system  may  be 
rendered  unstable  but  jumps  will  not  occur  unless  the  parameter  path  exits 
from  the  cusp  boundary. 

Maxwell  and  Delay  Rules 

An  interesting  application  of  catastrophe  theory  is  presented  in  Figure 
C-2.  Two  rules  of  behavior,  or  choices,  are  indicated,  which  assume  that 
behavior  may  in  some  fashion  perform  as  a system  which  seeks  a condition  of 
stable  or  minimum  potential.  If  a policy  is  at  point  PQ  on  an  issue  which  is 
polarized  around  policies  Pj  and  the  Delay  rule  invokes  a smooth  jump  to 

policy  Pj  and  the  Maxwell  rule  invokes  a discontinuous  jump  transition  to 
policy  P£.  Under  the  Delay  rule  the  jump  is  attended  by  a smooth,  non- 
dif f erentiable  change  in  the  state  variable  which  changes  the  policy  in  the 
direction  which  locally  maximizes  support.  Under  the  Maxwell  rule  the  jump  is 
accompanied  by  a discontinuous  change  in  the  state  variable  which  changes  the 
policy  to  where  the  support  is  globally  maximized. 

The  application  of  this  idea  lies  in  identifying  regions  of  instability 
in  system  state  variables,  and  being  able  to  anticipate  changes  in  the  state 
variables  by  designing  appropriate  system  response  functions  to  signal  system 
instability,  or  to  detect  that  a change  in  policy  is  required  to  maximize 
ef fectness . 
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Figure  C-l  : Catastrophe  Cusp  Geometry 
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